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ON THE SEPARATION OF GAS MIXTURES BY 
SUCTION OF THE THERMAL-DIFFUSION 
BOUNDARY LAYER 


By V. C. LIU 
(The University of Michigan, Ann Arbor, Michigan) 


[Received 21 November 1957] 


SUMMARY 


This paper discusses the formation and characteristics of the boundary layer of 
thermal diffusion existing along the surface of a heated wall over which flows a 
mixture of gases of unequal molecular weights. Due to the existence of a temperature 
gradient in the boundary layer, a thermal-diffusion flux is set up such that the lighter 
gas tends to move to the hotter region and the heavier gas to the colder region. 
When the temperature of the surface is maintained at a much higher level than 
that of the free stream, the concentration of the lighter gas increases monotonically 
from its free-stream value at the outer edge of the boundary layer to a maximum 
value at the plate. 

The existence of an extremely high temperature gradient makes it possible to 
obtain significant diffusive separation of gas mixtures in the boundary layer of 
a laminar flow. From this consideration a new method of separating gas mixtures 
by suction of the thermal-diffusion boundary layer was conceived. An analysis of 
the thermal-diffusion boundary layer is made, using a flat plate with constant 
suction as a model and assuming that the concentration of the lighter gas is much 
smaller than unity. The strong stabilizing effect of suction on the laminar flow is 
discussed briefly. A limiting ratio of suction against free-stream velocity for the 
maintenance of laminar flow is derived. An estimation of the rate of attaining the 
equilibrium-concentration profile is also made. 


1. Introduction 


Tue Clusius and Dickel method (1) for separating gas mixtures with an 
appropriate combination of thermal diffusion and natural convection has 
proved to be exceedingly effective in many cases. The principle of this 
method can be described as follows: a vertical tube has an electrically 
heated wire along its axis and is filled with a mixture of two gases of 
different molecular weights. On account of the temperature difference 
between the heated wire and the cooled wall, the resulting thermal- 
diffusion flux makes the concentration of the heavier constituent greater 
near the wall and that of the lighter constituent greater near the wire. 
Natural convection produces a flow of the cold gas downward near the wall 
and of the hot gas upward along the wire. These two processes combine 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 1, 1959) 
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2 Vv. Cc. LIU 
to increase the concentration of the heavier constituent at the bottom 
and of the lighter constituent at the top of the tube. 

From the hydrodynamic point of view, a limitation of this separation 
method by convection and thermal diffusion is that the critical Reynolds 
number of transition from laminar to turbulent flow must not be reached. 
According to Onsager and Watson (2), the magnitude of the transition 
Reynolds number of the Clusius and Dickel separation column is of order 
one-tenth the values obtainable for pure pressure flow through pipes. As 
a result of the low transition Reynolds number, the efficiency of produc- 
tion of concentrated material is limited (Jones and Furry, 3). 

The purpose of this paper is to point out a method of separating gas 
mixtures which might possibly be a more efficient method because of its 
inherently higher transition Reynolds number and higher rate of separa- 
tion. The method originates from the theory of laminar boundary layers, 
a phenomenon of fluid dynamics. It is most simply described by stating 
that if a forced convection stream at high Reynolds number flows along 
a stationary solid wall which is maintained at an elevated temperature, 
a very thin layer of laminar flow exists in contact with the wall. It is in 
this layer that the velocity of flow rises from zero at the wall to its free- 
stream value at the other edge of the boundary layer. Coexisting with 
the velocity boundary layer is a thermal boundary layer in which the 
temperature difference (7'—T,,) varies from (T,—T,,) at the wall to zero 
at the free stream. With a flow of ordinary gases, the thicknesses of these 
two boundary ‘ayers are of the same order of magnitude. For example, 
a temperature gradient of the order of 3 x 104 °C per cm can be obtained 
at the inner edge of the boundary layer with a forced convection of air 
under normal conditions, at a free-stream velocity of 30 m per sec, and 
a wall temperature of 300° C above the free-stream air temperature. 
Therefore, one might expect that if the flow medium is a mixture of gases 
of unequal molecular weights, diffusive separation will occur in the 
boundary layer, especially near the wall. The amount of separation after 
equilibrium is attained will depend on the local temperature gradient, the 
molecular-mass ratio, and the law of mutual interaction of the gases. 
Concentration of the lighter constituent will increase near the wall. This 
increase of concentration over the free-stream concentration diminishes 
to zero as the other edge of the boundary layer is approached. Thus, a 
concentration (or thermal-diffusion) boundary layer is formed. To collect 
the more concentrated mixture near the plate, one may use the method 
of boundary-layer suction through a porous plate, a technique already in 
use for the control of boundary-layer flow in attaining the optimum lift 
and drag characteristics of an airfoil (4). 
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The Clusius and Dickel method differs from the present one in the use 
of the convection stream and the collection of the concentrated mixture. 
The convective flow in the thermal-diffusion column of Clusius and Dickel 
is caused solely by the density gradients created by the temperature 
differences in the horizontal direction. The convective flow in the system 
here proposed is created by external causes (e.g. a pump). Such a flow 
is generally called a forced convection as distinguished from the natural 
convection caused by temperature differences and buoyancy forces. The 
forced convection is directly responsible for boosting the temperature 
gradient in the boundary layers. 

Another asset of the present method is its inherently higher transition 
Reynolds number due to the large stabilizing effect of suction on the 
boundary-layer flow. The stabilizing effect of suction makes it possible 
to use a higher temperature difference for thermal diffusion. 

The analysis that follows is to establish the validity and feasibility of 
the present method. The emphasis will be mainly on the concepts and 


on the general features of the problem rather than on the computation 
of details in specific cases. 


2. Notation 
The notation of Chapman and Cowling (5) has been used with the 
following exceptions: 
velocity vector of mean mass flow. 
’ mean mass-flow velocity components parallel to the x- and y-axes. 
kinematic viscosity (y/p). 
ordinary coefficient of diffusion of binary gas-mixtures. 
(7 —T..)/T x. 
logarithmic derivative of the gamma function. 
Prandtl number (yc,/A). [o is assumed constant for a given gas or 
gas mixture. ] 
v D (|S is assumed constant for a given gas mixture]. 


A/pc, D (A is assumed constant for a given gas mixture]. 
Mzy/(M,—M,). 


uv 
suction boundary-layer parameter E = —Py Up | (1/p) ay}. 
0 


viscosity- and diffusion-coefficient exponent 
p/po = (T/T): nD/(nD),. = (T/T). 


separation efficiency (see equation (5.1)). 
(1/T)\(@T /eZ). 
Mach number. 





¥. & 


R_ Reynolds number (uw, L/v..). 
L acharacteristic linear dimension of the system. 
R, u,,8*/v, (transition Reynolds number). 


Subscripts 
9 _—sc at: alll. 
at free stream. 
at suction chamber. 
constituent ‘1’ of a binary gas mixture. 
constituent ‘2’. 


3. Convective-diffusive equation 
Consider a non-uniform binary gas mixture. Let C, denote the mean 
molecular velocity of constituent 1 relative to a coordinate system 
moving with the mean mass flow velocity V. By applying the principle 
of conservation of molecules to the first constituent, one obtains 
on, 


+V.(n,V +n, €,) 0. (3.1) 


ot 


On adding (3.1) to the corresponding equation for the second constituent, 
and noting that p,€,+p,€, = 0, one gets the equation of mass conserva- 
tion for the mixture 


P4V.(pV) = 0. (3.2) 
ot 


For definiteness, suppose m, > m,. Let @ = n,/n and B = m,/(m,—m,). 
Substituting @ and 8 in (3.1) and (3.2) and then eliminating p, one has 


n = nv.VO+ PPV. (nC) = 0. (3.38) 
a 


The diffusive flux, given by 


“aS In p—af(1—6@)V in TI, (3.3 b) 


n,C, = —nD|V04 
includes contributions from gradients of concentration, pressure, and 
temperature, respectively; here a is the thermal diffusion factor. The 
diffusive flux due to external forces is taken to be negligible. The formula- 
tion of the total diffusive flux in (3.3) is based on the Chapman—Enskog 
theory of diffusion (5, note that the lighter gas is labelled 1 here). 

The diffusion coefficient D is inversely proportional to n if the tempera- 
ture is constant and nD is a function of T only. For many gas mixtures 
it can be assumed (5) that 


nD/(nD),. = (T/T..)*, (3.4) 


where s is a constant depending upon the molecular models employed. 
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For ordinary gases s lies between the extreme values } and 1. The value 
of s for a mixture of H, or He in air is 0-77. 

The dependence of the thermal-diffusion factor « on 7’ and n,/n, is 
small and here neglected. A list of measured values of « for various 
mixtures of gases has been compiled by Grew and Ibbs (6). Strictly, the 
functions V, n, and 7 appearing in the differential equation (3.3) for 4 
depend, in turn, on 6. This dependence is negligible when @ is small. 


4. Boundary layer with suction 

The boundary layer theory deals with the velocity and temperature 
distributions in the immediate neighbourhood of a solid surface along 
which there is convective flow of a slightly viscous medium. In the case 
of a flat plate with constant and continuous suction distributed over the 
whole plate, the Navier-Stokes equations of motion possess a particular 
solution for which the velocity distribution is independent of x, a coordi- 
nate measured along the plate (4). It can be shown for the case of an 
incompressible flow (7) that this special velocity profile, known as the 
asymptotic suction profile, is attained after an initial plate length of 


Xr, = 4u,,v/r? (4.1) 


measured from the leading edge, vy being the normal velocity at the plate. 
The corresponding problem for a compressible flow with heat transfer 
was treated by Liu (8). A brief résumé of the results of the asymptotic 
suction profile considering heat transfer and compressibility is reproduced 
here since the cited reference is unpublished. The equations of motion, 
continuity, and energy of boundary-layer flow along a flat plate are 


respectively (9) Cu ou é ou 
CP 4.2 
OM oa t PP ay a(t A il 


and sO) +5 (00) == 0, (4.3) 


3 @ é(,é éu\? 
7 é _ 6 (,eT\, (eu)? 4.4 
and a (cy T)+pv oy (7) ey ( ay) (ey) re 


The associated boundary conditions for the case of a flat plate with 
constant suction v, are 


-0, v=-—w, T= aty= 0, (4.5a) 


and e=u,, T=T at y = o. (4.5b) 


In view of the asymptotic assumption, i.e. @u/éx = 0 and @7'/éx = 0, one 
can integrate (4.3) immediately, giving 


pv = constant = po) = pu Va (4.6) 
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For many gases, the variation of » with 7 can be adequately represented 
(|b) = (T/T) (4.7) 
where the values of s lie between } and 1. For air s is about 0-76. 
The values of s in (3.4) and (4.7) are assumed identical. 


Introduce and 
substitute 


Po | (1/4) dy (4.8) 
0 


in (4.2) and (4.4). One obtains after simple calculation 


u u,(1—e-$), (4.9) 
T \ ae , . 
and —_ i} 7% +-4(y—1)M%, ~ (e- %—e-*s), (4.10) 
i. 2—o 
where the Prandtl number o = yc,/A is assumed constant. It should be 
noted that the corresponding solution in the case of an insulated plate 
was given by Young (10). 
It is significant that the essential variation of u and T' is restricted to 
a very thin boundary-layer zone, the thickness of which is of the order 
of v (—v»). Hence, it is justifiable to introduce the usual boundary-layer 
assumptions wu > v and @/éy > é/éx. For the case 6 < 1, the convective- 
diffusion equation (3.3) is simplified (by boundary layer approximations) 
to 


j } } ain T" 
.— { nv” a. inp (= 2m “)} (4.11) 


ct x cy oy | oy oy 


After the velocity and temperature distributions in the boundary layer 
attain their respective asymptotic profiles, (4.9) and (4.10), along a flat 
plate with constant suction, the concentration distribution @ will adjust 
itself suitably. 

Assume now that the origin of the x-axis coincides with the point at 
which the asymptotic suction profiles have been essentially attained (see 
(4.1)). Lett a u,.t; then by substituting (4.8) in (4.11), this becomes 

pus : ~\< paral 84 “(er 4... “)| _ 0. (4.12) 
u.,} Ct G o\og a 


where S = v/D, which is assumed constant for a given gas mixture. 


5. The steady-state solution 


Consider a Hat porous plate, on one side of which is a suction chamber. 
The suction-chamber pressure is maintained at a constant value p,, which 
is less than the free-stream pressure p, on the other side of the plate. 
A binary gas mixture, with @ < 1, flows along the plate with a free-stream 


+ This is feasible because the asymptotic profiles are independent of z. 
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velocity u.. The plate temperature is maintained at 7, which is higher 
than the free-stream temperature 7,,. It is further assumed that the 
composition of the gas in the suction chamber remains constant, i.e. 
6. = constant. 

The gas mixture, while seeping through a porous plate or an equivalent 
array of capillary tubes, is subjected to a diffusive-separation process. 
For a given p, and p,, the diffusive separation factor 

6.(1 —O9)/Og(1—O,) = 8,/05 
depends primarily on the ratio of the mean free path of the inlet mixture 
1, to the radius r, of the equivalent capillary tubes. When / >r, the 
flow through the porous plate obeys the well-known Knudsen formula 
(11); hence, @,/0, at p, = 0 has its maximum value (m,/m,)'. As the ratio 
l/r decreases, the intermolecular collisions ~vhich transfer momentum 
between the lighter faster and the heavier slower molecules become more 
frequent; consequently, both constituents of the mixture tend to move 
with a common mass velocity. As 1 <r, the flow pattern gradually con- 
forms to viscous tube flow of the Poiseuille type, for which the diffusive- 
separation effect is insignificant. In the intermediate region, 0,/0, varies 
between these extreme values. Let the separation efficiency Z be the ratio 
of the actual enrichment (increase in concentration) to the ideal enrich- 
ment at constant back composition, namely (12) 
0.—9% a EZ (m,/m,)*¥—1 , (5.1) 
9 (m,/m,)* 
where Z depends on po, p,, and l/r (see Appendix). 

The best porosity size for the experiment depends on the limiting suction 
velocity for stabilizing the boundary-layer flow (see Section 7), the separa- 
tion efficiency Z, and economy of operation. An optimum value for r 
can be determined from these considerations. The use of a dual system 
of suction, a combination of porous surface and slit orifices, might possibly 
provide the optimum condition for efficient separation. 

In the case of the steady state, (4.12) becomes 


26 dé g@in 
wet Saar ar) = ° 
with the boundary conditions 
6=0,; d0/di=0 atl =o. 
After elementary integration, (5.2) and (5.3) give 


6(£)/0,, = 1—aT'%e-8 f eS (dT dé) T 2+ dé. 
t 
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The diffusive separation factor for the thermal-diffusion boundary layer, 
6,/@,., in the case of incompressible flow (M,, < 1), becomes 


, 1 
5 1+ aa(1+-a)* [ 2 A(1+az)-4+® dz, (5.5 a) 
el) 0 

provided 0 < A < 1 where A = S/o which is constant for a given gas 
mixture. It is of interest to note that the definite integral in (5.5a) can 
be expressed in terms of hypergeometric functions (13); consequently, 


6 x 
b 


5° 1+-aa(1-+-a)*(1—A)-1F(1+a,1—A; 2—A; —a). (5.5b) 
The advantage of the form (5.5b) is that in special cases a simple formula 
for [@,/6,.|, can be obtained in closed form (in view of the considerable 
amount of knowledge concerning special forms of hypergeometric func- 
tions). 

The factor A, defined as A/pc,, D, needs further explanation. The para- 
meter A/pc,, which appears in Fourier’s equation of unsteady heat flow is 
sometimes called the thermal diffusivity (following Kelvin). It has the 
dimension L?/t which is the same as that of the mutual-diffusion coefficient 
D; hence, A is dimensionless. Since both A/pc, and D are proportional to 
T+* (5), A is nearly temperature-independent. 

Values of A can be estimated on the basis of the kinetic theory of gases. 
For isotope mixtures, the mutual-diffusion coefficient must be nearly 
equal to the coefficient of self-diffusion of either isotope. Using the self- 
diffusion formula for Maxwellian molecules and Eucken’s formula for 
thermal conductivity, one obtains A = 0-9 for diatomic mixtures. The 
value of A for a mixture of H, and N, is estimated to be 0-32, using the 
mutual-diffusion coefficient. 

In the case of an isotope mixture, « is usually very small compared with 
unity. If a = 1, a simple formula for the diffusive separation factor for 
the boundary layer can be obtained, namely 


& 1+-a2%-My(1—$A)—J(}—$A)]- (5.6) 


“b 


On the basis of (5.6), one can calculate the enrichment factor (@,—9@,,)/@... 
This is found to be about fourteen times that corresponding to the con- 
vection-free case, namely (7,/7,,)*—1, when the value of « lies in the 
range 0-01 to 0-1. 

In view of (5.1), one obtains the total diffusive separation factor 


fee} = [4G fet), 
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6. The time to attain the steady state 

The time required to attain the steady state will largely determine 
whether or not this method of separating gas mixtures by suction of the 
thermal-diffusion boundary layer is convenient and practical. It will be 
mathematically difficult to determine the time lag of the system by 
solving the time-dependent equation (4.12). A careful investigation of the 
problem, however, reveals that most of the concentration change occurs 
near the wall surface of the boundary layer; hence it u:ay suffice to solve 
an approximate form of (4.12), valid for the region near the plate. The 
temperature distribution near the wall surface of the boundary layer can 


be approximated by T “= Ty exp(G, £), (6.1) 
where G, is determined by matching the temperature gradient of (6.1) with 
that of the exact expression (4.10) at ¢ = 0; hence 

Gy = [dT /dl),/T). (6.2) 


The velocity near the plate is taken as u = 0. With these approxima- 
tions, (4.12) becomes 


00/ét = (v3 /v_ S)[6*0/el?—(aG,—S)00/0C}, (6.3) 


provided one assumes that yz/u, = 7'/T). The boundary conditions are 


6(0,Z) = 0,, A(t, 0) = ,, 0(t,5) = 05 = 6... (6.4) 
The simple transformation 
6 = O* exp[4(aG,— S)l—(aG,— S)*v§ t/4v,. 8] (6.5) 
applied to (6.3) leads to 
20* /et = (v3/v, S)(620*/at?), (6.6) 
a form of the familiar heat equation. The transformed boundary condi- 
tions (6.4) are as follows: 
6*(0, ¢) — 6. exp[—4(aG,— S)¢], (6.7 a) 
0*(t,0) = 0, exp[(aG,— S)*v3 t/4v, 8], (6.7 b) 
6*(t,5) = Pexp[— }(aG,— S)5+ (aG,— S)*v§ t/4v, 8]. (6.7 ¢) 
The solution to the boundary-value problem of (6.6) with (6.7) can be 
constructed by the classical method of Laplace transformation or Duha- 
mel’s theorem. One may note, however, that the present boundary-value 
problem is completely analogous to that of heat conduction in a solid bar 
with initial and boundary conditions prescribed as in (6.7) (14). 
Using this mathematical analogy, the concentration near the plate is 


found to be 
A(t, f) = O* exp[4(aG,— S)f—(aG,— S)*v3 t/4v, S], 
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where 


4n75 sin(nnl/8) 
1 5*(aGiy— 8)? dnt? 


O*(t, C) (2/8) ((0..- -Oy)e-M*ttv*tiveSd* | 


L {8,—(— 1)"6, € (aGe—S)8/2} o(aGe—S)rgtiaveS) (6.8) 

To determine whether the present system will take an unduly long time 

to attain the steady state, it is necessary to estimate the relaxation time 

of the system by taking merely the first eigenvalue n = 1 in (6.8). This 
gives a relaxation time 

t, ~ 4v, S8?/v2[ 40? +-8?(aG,—S)?]. (6.9) 

It is apparent from (6.9) that with values of S and 8 both of order unity, 

the governing factor for t, is vg/7?vz. This is of order 0-1 sec if vy is 1 cm/sec. 


7. The stabilization effect of suction on the boundary layer flow 

The present analysis is supposedly valid only in the case of a non- 
turbulent boundary-layer flow, although in the case of turbulent flow a 
laminar sub-layer along the plate still exists. It is generally known, how- 
ever, that turbulent motion is the more natural state of fluid motion, 
whereas laminar motion occurs only when the Reynolds number is so low 
that deviations from it are liable to be damped out. The transition 
Reynolds number for an incompressible boundary-layer flow depends upon 
such factors as the free-stream turbulence, the surface temperature, 
pressure gradient, and the roughness of the surface, to mention a few 
(15). 

Suction reduces the boundary-layer thickness. A thinner boundary 
layer is less prone to become turbulent; furthermore, it creates a laminar 
velocity profile which has a higher limit of stability than that without 
boundary suction. 

Experimental investigations with stability analysis (16, 17) appear to 
confirm the suggestion that above a certain value of —v,/u,, against in- 
finitesimal disturbances, complete stability of the isothermal laminar 
boundary layer along a flat plate with constant suction may be achieved 
for all Reynolds numbers. Liepmann and Fila (18) provided experimental 
results of surface temperature on the boundary-layer transition without 
suction. Although precise information concerning the stabilizing effect 
by suction on the boundary layer along a heated plate is still lacking, 
from the expression for the displacement thickness of the boundary layer 
along such a plate 


@ 





ON THE SEPARATION OF GAS MIXTURES ll 


one ean glean a suggestion of the stabilizing effect exerted by suction. 
(This, of course, does not disregard the known fact that, in some cases, 
surface temperature can have more important effects} on the stabilization 
of flow.) From (7.1) and the definition of R,, one obtains 


x See fe T..\'~* 7.2 
vo = ot { (a )FR) (7.2) 
0 


After substituting the asymptotic suction profiles (4.9) and (4.10), this 
can be integrated numerically. In the special case s = 1, (7.2) can be 
integrated in closed form, 


— 1 _[1(To_)\, 2-1) 1 7 
= BT. “(7 +500 2m," — 3) +1]: (7.3) 


In the case of isothermal incompressible flow, the suction ratio is 


—v,/u.. = 1/R, (7.4) 


which is of order 10-4. 

It should be noted that the asymptotic suction profiles, on the existence 
of which the above calculations are based, develop only after a certain 
distance from the leading edge (see (4.1)). The velocity and temperature 
profiles for the intermediate region have lower limits of stability, hence 
the suction ratio over this region must be larger than (7.3) if laminar flow 


is to be maintained. 


8. Conclusion 


The case of a flat plate with constant suction and @ <1 is chosen as 
a model for the purpose of analysis. The solution, simple as it is, pre- 
serves the general features of a thermal-diffusion boundary layer with 
suction. This solution, being exact, serves well to explore the mechanism 
of the physical phenomenon which is the primary aim of the present paper. 

The relaxation of such restrictions as 6 <1, dp/éx = 0, and % = 
constant would involve more laborious, but straightforward, calculations 
(19). 

The radiative heat transfer between the surface and the flowing gases 
has not been taken into account in the analysis. This would tend to make 
the present results of diffusive separation appear too high; but the excess 
is believed to be small unless the surface temperature is much higher than 
three times the free-stream temperature. 

Although the primary function of this discussion is to establish the 


+ These include the effect of gravitational forces (tecause of the density differences in 
the boundary layer) and the effect due to the dependence of the viscosity of a gas on the 
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principle of separating gas mixtures by suction of the thermal-diffusion 
boundary layer, it is of interest to suggest a simple scheme of application 
of this principle. To produce the forced convection of the mixture for 
separation, one may use a conventional low-turbulence pressure channel 
(20). Unlike the wind tunnel usually designed for aerodynamic tests of 
models, the pressure channel in question should have at its test section 
properly-spaced porous plates designed for the production and suction of 
the thermal-diffusion boundary layers. If the channel is geometrically 
similar to the Langley two-dimensional low-turbulence pressure channel, 
and can, for example, displace 5 kg of gas mixture per second at a speed 
of 30 m per sec at its test section, the power consumed in the channel will 
be about ? kW. This does not include the power lost in the drive-unit, 
nor that used in heating the plates, sucking the boundary layers, and 
cooling the mixtures in the channel. The power consumption is estimated 
on the basis of a power factor of 0-32, an average value for tunnels of this 
type. The above estimation is meant to give the right order of magnitude 
only. In order to achieve a large change in gas abundance, one needs to 
use a cascade of these processes. 
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APPENDIX 
CAPILLARY THEORY OF SEPARATION 


In the special case 6 << 1, the separation efficiency becomes (12) 
fo 
Z (Eo Ee, E) = | exp(€*) d&/(Eo—E-)expl&), 
where & = (128.S8/92r)#{1 + (92r/128)(r/1)] + (92r/512)4(.Sm,/m,)-4, 


E = (128S/97)§+- (97 512)4(Sm,/m,)-?. 
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HEAT TRANSFER BY LAMINAR FLOW FROM A 
ROTATING DISK AT LARGE PRANDTL NUMBERS 
By D. R. DAVIES (University of Sheffield) 

(Received 18 December 1957] 


SUMMARY 


In this paper the approximate method of calculating the distribution of rate of 
heat transfer by laminar flow from a flat plate, described recently by Davies and 
Bourne, is extended to the problem of a heated rotating disk for an arbitrary radial 
distribution of surface temperature on the disk. The analysis is applicable if the 
temperature boundary layer is embedded well within the velocity boundary layer 
over the disk. By comparing the ensuing numerical results for heat transfer with 
those calculated from an exact solution, by Millsaps and Pohlhausen in the special 
case of constant surface temperature, we find that very good accuracy is obtained 
when the Prandtl number is greater than about 6. 


1. Introduction 
MULLSAPS AND POHLHAUSEN (1) have recently obtained an exact solution 
(including aerodynamic heating) of the equation of forced convection from 
a rotating plate, when the flow is entirely laminar. This solution was 
used to compute the distribution of temperature and rate of heat transfer 
for values of the Prandtl number oc, between 1 and 10. However, the 
analysis is only applicable when the surface temperature is independent 
of radial and angular position on the plate. In order to obtain a solution 
when the surface temperature of the plate is an arbitrary function of radial 
distance, Morgan and Warner (2) have applied in a modified form the 
approximate technique used first by Lighthill (3) in the flat plate case. 
A linear representation of the velocity profile (valid in the region very 
near the surface of the disk) was assumed to apply over the whole thick- 
ness of the boundary layer. At very large values of o (o ~ 10 for cold 
water, o > 100 for oils) the thickness of the temperature boundary layer 
is very much smaller than that of the velocity boundary layer, so that 
the crucial part of the velocity profile for heat flow is that near the surface. 
Consequently a linear representation is likely to be adequate for large 
values of o. By comparing their numerical results with the exact values 
given by Millsaps and Pohlhausen, Morgan and Warner find that the error 
is large at o = 10, the highest Prandtl number for which heat transfer 
rates were computed by Millsaps and Pohlhausen, and they suggest that 
the smallest value of o for which the linear approximation is suitable is 
very considerably greater than 10. 

Davies and Bourne (4) have shown that this form of approximation 

(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 1, 1959] 
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can be improved by using a power law representation of the velocity 
profile and so produced a better approximation over a larger part of the 
boundary layer thickness than was possible by using the linear (tangent) 
representation. An application to the disk problem of the power law 
approximation is therefore likely to yield a very substantial reduction in 
the minimum value of ¢ at which good accuracy is obtained. 

In the problem of the heated disk the radial component of velocity 
conveys the heat in the air stream outwards and consequently plays the 
dominant role in the convection of heat. A very good power law repre- 
sentation of this component can be made for the inner one-sixth of the 
velocity boundary layer. Following Davies and Bourne (4), the basic 
temperature equation is considerably simplified by using a von Mises 
transformation with the power law representation. The method of sources 
is then applied, the heated disk being regarded as an assembly of con- 
centric circular sources. The radial distribution of source strength (or 
rate of heat transfer) is determined by solving an integral equation, the 
solution depending on the prescribed radial distribution of surface tem- 
perature on the disk. In the particular case of constant disk temperature 
the numerical values of rate of heat transfer for values of o greater than 
about 6 are found to be in very good agreement with the exact results 
given by Millsaps and Pohlhausen; these numerical values are given and 
discussed in section 4 of this paper. 

The advantage of this method of calculation lies in the comparative 
simplicity of the numerical evaluation and particularly in its application 
to cases in which the distribution of disk temperature is an arbitrary 
function of radial distance from the centre of the disk. An additional 
point of interest lies in the possibility of extending the analysis (as in the 
flat plate case, 4) to include the important practical problem of convective 
heat transfer by turbulent flow from a rotating heated disk. 


2. The basic temperature equation 
If r denotes radial distance from the centre of the disk, z normal distance 
from the surface of the disk, u, the component velocity in the radial 
direction, u, the component velocity perpendicular to the plate, v the 
kinematic viscosity, « the thermal diffusivity, 7’ the temperature, the 
equation of continuity is 
1 @ 
rér 
and the temperature equation is 


oT oT 1 é/2@ a 
“G+ oe = 5(S)+ Sl: 


ra 
(ru,)+ = %, = 0, (1) 
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The aerodynamic heating term has been neglected in equation (2) but this 
effect is independent of the distribution of surface temperature; it is 
dependent only on the velocity distribution and (as shown by Millsaps 
and Pohlhausen) can be added to the solution of equation (2). 


Az 





Fic. 1 


The distributions of wu, and u, are required in order to obtain solutions 
of (2) and these have been computed by Cochran (5). The tangential 
velocity component does not enter the temperature equation, since the 
problem is axially symmetric and tangential derivatives of temperature 
vanish. 

Following the usual boundary layer approximations, we neglect, on the 
right-hand side of (2), the terms involving derivatives of 7’ with respect 
to r in comparison with those involving derivatives with respect to z. 
This is particularly applicable in the case of the disk, as von Karman (6) 
has shown that the boundary layer thickness in the laminar case is 
approximately independent of radial distance. Millsaps and Pohlhausen 
have also shown, by calculation of the exact solution, that in practical 
cases this approximation is permissible. 

The basic temperature equation is then 

u or + U eo? = eet 


rs: Z 
cr 


(3) 


ez oz? 


This is further simplified by a von Mises transformation (see, for example, 
7, p. 126), taking r and & (the stream function) as independent variables 
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instead of r and z and using the continuity equation (1). We obtain 


é é 
— = «r?—ly, —]. 4 
a a a) 
The radial velocity u, has been computed by Cochran (5), and is written 
in the form u, = wrF(f), (5) 


where w is the angular velocity of the disk and { = (w/v)!z; a part of the 
computed function F(¢) is shown in Fig. 2. 
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Fic. 2. Caleulated values of F(f) = u,/(rw), where £ = (c/v)*z: 
— exact solution ; —-—-—- power law (F = 0-245{[°*”) ; 

— — — — linear law (F = 0-51{}*). 


It has not been found possible to apply the method of sources, used 
previously by Davies and Bourne (4) in the flat plate case, with the exact 
form of the function F. However, in many practical problems a is large, 
the temperature boundary layer is embedded within the velocity boundary 
layer and the method used by Davies and Bourne can be applied by 
choosing a power law representation in the form 

F = Am (6) 
to fit the u, profile in the inner range of { (ie. {<< 07). A very good 
representation is found to be possible (4 per cent deviation at most from 
the exact values in this range of {), and the ensuing calculated heat transfer 


5092.45 Cc 





18 D. R. DAVIES 


values are found not to depend critically on the choice of A and aq; i.e. 
any pair of values of A and a which yields a good fit to the actual profile 
for ¢ < 0-7 leads to calculated heat transfer values in reasonable agree- 
ment with each other. 
Using power law (6) and the continuity relation expressed in terms of y/, 
we obtain ab 
ae om FU, oi Aw) +t%y-tayp2z0 (7) 
and hence b = Awl+tdy-t( 1 + q)-tr2zt+a), 


taking y = 0 along z = 0. Expressing u, in terms of ys, the basic tempera- 


ture equation (4) becomes 
éT é/,,@ 
— = C— — ]|, (9 
a a" a) 


where R = rr2+a+a) 


t = a/(1+a), 


and c= hicey(* +2021 +a), boe(1+a) 4 UL +a)( ] 4 (1 +2ay(1 +a) 7 (2 1 x)~!, 


3. Solution for an arbitrary radial distribution of disk temperature 
Equation (9) is identical in form with the basic equation discussed 
previously by Davies and Bourne (4, equation 11). Suppose now that 
a continuous uniform circular source, radius ry, concentric with the disk 
and lying on the surface of the disk, emits Q units of heat per unit time. 
The analysis of the flat plate solution (4) can then be applied and a solution, 
giving the distribution of local rate of heat transfer, obtained for any 
prescribed radial distribution of surface temperature. If 7;(r,) and 7, 
denote respectively the temperature of the disk and of the stationary air, 
then the distribution of Q(r,), the local rate of heat transfer, is given by 
Q(r4) = m7? 2-28(1 + )?4(2 + a) 28-DDP(2/8)sin(22/s)A9-2k x 
ro 
aby bouasay S [ (rg— R)®*-2(T(r.)—T,] dR, (10) 
i eg 
where s = 2(2+-a)/(1+ a) and k is the thermal conductivity. 
This expression necessitates, in general, numerical integration and 
numerical differentiation, but in many cases the integral in (10) can be 
expressed in terms of tabulated functions or simple quadratures; e.g. if 


T,—T, = > a, 1", the integral can be written, as in the flat plate case 


n= 


discussed by Davies and Bourne (4), as a series of beta functions and the 
temperature distribution expressed in terms of a Pearson J function (8) 
and a series of quadratures (these have been shown recently by D. E. 
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Bourne, in a private communication, to be expressible as Whittaker func- 
tions). Once the distribution of Q(r,) has been evaluated, the distribution 
of temperature can of course also be evaluated, as in the flat plate case. 


4. Solution for constant disk temperature 


In the particular case of constant temperature, 7,, on the surface of the 
disk, the expression (10) for Q reduces to 


Q(ry) = 7-2 2-2(1 +-x)*(2 + x)**-1T(2/8)sin (27/8) A! kewtytg2+® 
x (T,—T,) r2-", 


The total rate of heat transfer from an area of radius r is given by 


E =| 27Q(ro)drp, 
0 
and we find on integration that 


E(k(T,—Ty)wtv-tr?}-? = Ko@+, (11) 
where K = 20-2024 )%*(1 + a)/8-DA4G-2/) (2 /8)sin(27/s); , 


we note that the form of the dependence of £ on r, w, and v is independent 
of «. The corresponding result given by Millsaps and Pohlhausen (1, 
equation 2.29) is 

E[k(T, —Ty)w'v-tr?]-? = 2Q}(0), (12) 
where Q)(0) is a function of o (1, Fig. 9). It is of interest to note here 
that for very large values of o, the temperature boundary layer is only 
a very small fraction of the velocity boundary layer; the (tangent) linear 
approximation to the velocity profile is likely to be adequate, we can then 
take « = 1, and equation (11) shows that the dependence of F on is given 
by o! as in the flat plate case. We note also that the asymptotic form of 
the Millsaps-Pohlhausen equation (12), as o-> 00, can easily be shown, 
using the first term of the series expansion given by Cochran (5) for u, 
and equation (2.19) of Millsaps and Pohlhausen (1), to coincide with 
formula (11), with a = 1 and A = 0-51; it is given by 

Elk(T,—T)wby-r2}-? = (0-51)'34T(§)sin(Gr)o! 

The numerical results calculated from formulae (11) and (12) are shown 
in Table | for various Prandtl numbers. 

The numerical values of rate of heat transfer given in the second column 
of Table 1 were computed from the exact solution (12), for o > 10, using 
the Millsaps—Pohlhausen expression for Q}(0) (1, equation 2.19), with the 
series and numerical values for u, given by Cochran; in this way the 
numerical values for EZ, evaluated by Millsaps and Pohlhausen for 
1 < o < 10, are extended from o = 10 up to o = 1,000. These numerical 
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results show that the approximate formula (11) with a power law repre- 
sentation (a = 0-67, A = 0-245) can be applied with very good accuracy 
for values of o between 6 and 200 (this range will include many cases of 
practical importance). For values of o less than 6 an effective part of the 
temperature boundary layer clearly lies outside the inner zone 0 < ¢ < 0-7. 


TABLE | 
Calculated values of rate of heat transfer from a rotating disk using the 
exact and approximate solutions 





Percentage over-estimates of 
Calculated values of E{k(T,—T,)r?wtv-*}™ the exact values by the approxt- 
in cm. sec. units mate formulae 





Using Using 
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At o = 200 we find that the asymptotic solution leads to the same error 
(4 per cent) as the approximate power law formula (11). If o > 200, 
however, the error due to the power law formula increases, and the error 
involved in applying the asymptotic solution of course decreases. The 
error due to the application of the linear approximation is considerably 
greater than that due to the power law approximation if ¢ is less than 200, 
but at higher values of o the linear approximation is more accurate than 
the power law approximation. However, an application in this range of the 
simpler asymptotic solution leads to an equally small percentage error. 
These results show, that we can apply the power law formula (11) with 
good accuracy for values of o between 6 and 200, while the simple 
asymptotic formula is sufficiently accurate for o > 200. 

A convenient expression for the temperature distribution is also easily 
derived, as in the flat plate case. We obtain 


= 1—I[2A(1+a)-(2+a)ef?+™, —(14+a)(2+a)-1], (13) 
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where 7 is a Pearson function (8). We again find that, for values of o 
greater than about 6, the temperature distributions calculated from (13) 
are in agreement with the exact values calculated by Millsaps and Pohl- 
hausen (1, Fig. 7) but the deviations increase as o decreases below o = 6, 
as in the calculation of heat transfer rates. 

We can, therefore, reasonably infer that for values of o greater than 6 
this approximate method of calculation is likely to be sufficiently accurate 
to calculate rates of heat transfer and temperature distributions for an 
arbitrary radial distribution of disk temperature, the associated numerical 
work in general consisting of the numerical integration and differentiation 
involved in expression (10). 
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VIBRATIONS OF BEAMS 


Ill. SCREW MODES+ 


By W. A. GREEN (Atomic Weapons Research Establishment, 
Aldermaston, Berks.)t 


[Received 14 November 1957] 


SUMMARY 
The exact equations of wave propagation in an elastic medium are solved to 
obtain the dispersion equation for the screw form of vibration of a circular cylinder. 
Dispersion curves are obtained for the fundamental mode, for a range of values of 
the Poisson’s ratio, in the long and intermediate wavelength region and the cross- 
overs of these curves and the corresponding curves for the longitudinal form of 
motion determined. 
Notation 
displacements in the directions r, @, z respectively. 
functions of r, see (1). 
arbitrary constants. 
Lamé’s constants. 


Poisson’s ratio. 


velocity of dilatation waves in an infinite medium. 


velocity of shear waves in an infinite medium. 
velocity of Rayleigh surface waves. 

phase velocity of waves in the bar. 

angular frequency. 

wave number. 


radius of bar. 


1+-A(z?+-y?)/2uy?. 


+ Part of a thesis submitted to the University of Wales for the degree of Ph.D. 
t Now at M.E.R.L., East Kilbride, Lanarks. 
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J2(x) 
Se 
J,(x) 
Jy) 
6, = y=. 
2= IT iy) 
6,(s) = 2—}s? (s <1). 
6,(8) = scot|s— }z], s real 
—is—}, s imaginary 


See Hudson (11). 


| \a| >. 


1. Introduction 

Ix the previous papers (1 and 2) the perturbation theory developed to 
determine the dispersion curves for elastic wave propagation in beams of 
non-circular section was found to break down at certain points in the 
intermediate wavelength region due to degeneracy between the screw form 
of vibration (3) of a circular cylinder and, on the one hand, the longitudinal 
form (1), on the other, the torsional form (2). This cross-over of the 
dispersion curves for beams possessing fourfold symmetry has been inter- 
preted as implying that in this region the screw becomes the fundamental 
mode of motion (4), and a further examination of this form of motion was 
felt to be worth while. 

The simplest way of examining wave propagation in elastic beams using 
the exact equations of motion is by considering beams of circular section 
where the boundary conditions are such that the different possible forms 
of motion can arise independently. In this way the longitudinal forms 
of motion have been considered in the intermediate wavelength region by 
Bancroft (5), Davies (6), and others (see 7), and in the short wavelength 
region by Adem (8), McSkimin (9), and Redwood and Lamb (10). In 
particular Bancroft has obtained the dispersion curves for ’ range of 
values of Poisson’s ratio for the fundamental mode of propagation in a 
circular cylinder, Davies obtaining similar solutions for the fundamental 
mode and first two harmonics for a Poisson’s ratio o = 0-29. The flexural 
form of motion has been examined in a similar way by Hudson (11) and 
Abramson (12), the former also obtaining some information about the 
higher forms of motion in the short wavelength region. In the following 
the dispersion equation for screw vibrations of a circular cylinder is solved 
numerically to obtain dispersion curves for the fundamental mode over 
the intermediate wavelength range for a set of values of the Poisson’s 
ratio analogous to the solutions of Bancroft and Hudson for the longi- 
tudinal and flexural vibrations respectively. 


2. Dispersion equation 
The Pochhammer-—Chree equations governing the propagation of elastic 
waves in cylindrical coordinates have solutions, for waves travelling along 
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the axis of z, of the form 


u = u,,(A,,cosn6+ B, sin n0)cos(kz—wt) 
v = v,,(A,,sinnO— B,, cos né)cos(kz—wt) (i = 1, 2,3) (2.1) 


w = w,,{A,,cosnb+ B, sin n@)sin(kz—wt) 


where the w,,;, Ung, Wp, are functions of r only and are given elsewhere (1) 
and the A,,, B,,; are arbitrary constants. The general solution for the 
displacements consists of linear combinations of these solutions over all 
values of n, but as has previously been pointed out (1 and 11) the boundary 
conditions for free vibrations of an infinitely long circular cylinder are 
satisfied by displacements involving each value of n separately. In fact 
for any given value of n the boundary conditions are satisfied by displace- 
ments involving only the A, or only the B,,;, one set being obtained from 
the other for all n > 0 by a rotation of axes through an angle 7/n. 

The screw form of motion is characterized by two nodal planes through 
the axis of the cylinder, the corresponding displacements being given by 


(2.1) with m = 2. Imposing the boundary conditions leads to the disper- 
sion or frequency equation 


{xJ3(x)—J,(x)}/x? $J3 (x) i{J3(y)—(K— 1)Ja(y)} 0 
{xJ a(x) — (hx? — 4)J,(x)}/ 22? { Jo(x)—2xJ5(x)}/x? {J(y)—yJa(y)}/y* 
—kaJ,(x)/2x* (k?— 8?) J3(x)/4kB —kJ3(y)/2a 
(2.2) 
or in terms of the @, functions defined by Hudson 
| 2(@.—1) 4—2°—0, 4—Ky*—6, | = 0. 
6,—4+}22 2(1—6,) 2(1—6,) | 2.3) 


2 (1—a?/z?)0, 26, 


3. Dispersion curves 


For any given value of the Poisson’s ratio o equation (2.3) may be 
solved to give the dispersion curve relating the reduced phase velocity 
c/c, to the number of wavelengths in the circumference of the bar ka. 
An infinite number of such curves exist for any value of o, the solution 
giving the smallest value of the phase velocity for any ka being termed 
the fundamental mode of motion, that giving the next value the first 
harmonic and so on. We consider the solution for the fundamental mode 
of motion only. 

Assuming a solution for which the phase velocity remains finite and 
bounded as the wavelength is increased indefinitely, both z and y become 
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small and with the appropriate expansion for the @, function equation 
2.3) may be written 


Pt 

aly") = 0 (3.1) 
correct to second order. This has solutions 
or r=Yy, ¢ = 9. 


Of these, the first satisfies equation (2.3) identically and for this case the 
displacements are no longer given by (2.1), whilst the second solution is 
trivial. Thus no solution of the screw dispersion equation exists for which 
the phase velocity remains finite and bounded at large wavelengths, and 
the velocity of propagation of the fundamental mode must therefore 
become infinite as ka tends to zeru, 

In the intermediate wavelength region the dispersion equation is simpli- 
fied to a greater or lesser extent at the points c/c, = 1 (x = 0); c/eg = v2 
(x = z), and c = ¢, (y = 0), and solutions were accordingly obtained at 
these points first, for a Poisson’s ratio o = 0-3. This allowed a sketch 
curve to be drawn which served to provide initial trial values for the 
numerical solution of equation (2.3) at any value of ka, the solution then 
proceeding by iteration. The exact curve for o = 0-3 having been ob- 
tained this in turn provides the initial trial values for other values of oc. 

The solution in the short wavelength region has been investigated for 
the general dispersion equation (n) by Hudson (11), who shows that the 
limiting value for velocities below the shear wave velocity c, is that of 
the Rayleigh surface waves c, in the medium and this has been verified 
for equation (2.3). 

The results for the fundamental mode of screw vibration are given in 
Table 1 for six values of the Poisson’s ratio over the range ka = 1 to 5, 
and the dispersion curve for o = 0-3 plotted in Fig. 1 together with the 
corresponding curve for the longitudinal form of motion obtained from 
Bancroft’s results.. Of particular interest are the cross-overs between the 
screw and longitudinal curves in the intermediate wavelength region, and 
these have been obtained graphically using Bancroft’s solutions for the 
longitudinal curves, extrapolated to o = 0 and o = 0-5, and are given in 
Table 2. 

That these solutions do in fact represent the fundamental screw mode 
has been verified by examining the cross-over of the curves and the line 
c = ¢,. At the shear wave velocity (z = 0) equation (2.3) reduces to 

24+ 1227+ 144 


6, = eat) (3.2) 
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Fic. 1. Dispersion curves for (a) longitudinal, (b) screw, modes (¢ = 0-3). 


TABLE 1 


Phase velocity as a function of wavelength and Poisson’s ratio 
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TABLE 2 


Intersections of longitudinal and screw dispersion curves 





or2 0-3 





1-680 1-714 | 4°645 
2°100 | 0-967 
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Considering the equation 
fie) = 8y)--- a 
since y is imaginary, f(z) is a regular function of z with the limits 
f(z) = (1—e/e})z*/6+-O(z*) (z+ 0), 
f(2) = —2*/6+ O(z) (z > 20), 
and there exists at least one real positive root of equation (3.3). Denoting 
the least real positive root of (3.3) by z, we must have that 
df —4—8+2i(1—cF/ef) (28+ 24) _ 
dz, Zy 3(22+-12)? ~ 
which on substituting for 6 from equation (3.3) leads to the condition 
z+ 3624 c3/c? — 288z?( 1 — 3c3/c?) —5184(1—c3/c?) > 0. (3.4) 





0 


This is satisfied provided z, > z’ where z’ is the only real positive root of 
the equation in (3.4). Further if the equality holds then we have 


Gf  xA(z{+36z2{+ 2882}+ 6912) (3.5) 

dst 18(z?+-12)8 ; 
and z, = 2’ is a maximum of the function f(z). This is clearly impossible 
and hence we must have z, > 2’. If z, is any other real positive root 





of equation (3.3) then since z, >z, we would again have df/dz, < 0 
showing that no such root z, exists. Thus there exists one and only one 
dispersion curve for any given Poisson’s ratio which crosses the line 
¢ = ¢,. Since no solution of the dispersion equation exists for which the 
phase velocity is always less than the shear wave velocity then the curve 
which crosses the line c = c, must correspond to the fundamental mode 
of motion, and this is the solution which we have picked out. 


4. Discussion 

The results given in Table 1 may be interpolated to give the dispersion 
curve for the fundamental mode of screw vibration of a circular cylinder 
of a material of any Poisson’s ratio in the long and intermediate wave- 
length region. In the region of short wavelengths and for higher modes 
an approximate method (13) has been employed to give qualitative and 
partly quantitative results. 

From the solutions of the determinantal equation (2.3) it is possible to 
obtain the ratios of the A,, or B,, and thus determine the displacements 
corresponding to this mode of motion. This has not been carried out here 
but the general behaviour of the displacements may be inferred from 
equations (2.1) and have been sketched in Fig. 2 over two cross-sections, 
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half a wavelength apart. This brings out the squeezing or screw character 
of the motion and the fourfold symmetry involved. It may be seen that 
the two axial planes which are nodal to the displacements in the longi- 
tudinal and radial directions are antinodal to the transverse displacement 
and vice versa. Displacements of this type would in practice be difficult 
to excite by the usual quartz crystal driving methods but a suitable means 
of squeezing together the opposite faces of a square bar might be expected 
to give rise to this form of vibration. 
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Fic. 2. Screw mode displacements for circular cylinder over two cross-sections 
half a wavelength apart. 
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SUMMARY 


By means of the Rayleigh—Ritz method, an analysis is made of the vibration of 
a rectangular plate whose edges are elastically restrained against rotation. Plate 
defiexions are represented by a set of functions which define the normal modes of 
vibration of a beam whose ends are elastically restrained against rotation. Values 
of various integrals of these functions and their derivatives are established. Fre- 
quencies are obtained from a set of linear simultaneous equations which may be 
solved by a simple iterative process. An approximate frequency equation is also 
derived and numerical tables for use with this equation are presented. 


1. Introduction 


Exact solutions to the differential equation of transverse vibration of a 
thin rectangular plate can be obtained only if one pair of opposite edges 
are simply supported (1, 2). For other boundary conditions approximate 
methods of solution must be adopted. Keference should be made to 
Warburton (3) for a comprehensive summary of these methods, as applied 
to plate vibrations, and to Hearmon (4) who quotes the numerical results 
of various authorities. However, all solutions, exact and approximate, 
known to the author, deal with plates under either simply supported, free 
or clamped boundary conditions, or mixed boundary conditions involving 
these three types of support. 

In the present paper we consider the vibration of a plate whose edges 
are rigidly supported against displacement, but which are elastically 
restrained against rotation. Opposite pairs of edges are taken to be re- 
strained to the same degree, although the method can be extended to a 
plate whose four edges are restrained to differing degrees. The Rayleigh- 
Ritz energy method is employed, and so the calculated frequencies will 
be somewhat higher than the exact values. Frequencies and modes of 
vibration are obtained from a set of linear simultaneous equations, which 
may be solved by a simple iterative process. An approximate frequency 
equation is also given, which is considered sufficiently accurate for 

[Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 1, 1959) 





30 T. E. CARMICHAEL 


practical purposes, and numerical tables for use with this equation are 
presented. 


2. Notation 
a,b lengths of sides of plate (b/a < 1). 
h thickness of plate. 
E Young’s modulus. 
v Poisson’s ratio. 
density of plate material. 
bending stiffness of plate Eh3/12(1—v*). 
gravitational acceleration. 
angular frequency, radians/sec. 
rectangular coordinate axes. 
amplitude of plate vibration. 
maximum elastic strain energy of plate. 
maximum kinetic energy of plate. 
boundary restraint parameters. 
span of beam. 
coefficient in series representation of w(x, y). 
characteristic beam functions defined by equation (7). 
constants of integration associated with F.. 
characteristic value of F.. 
a function of 2 alone. 
a function of y alone. 
Lim Nims Him) 
Ems Mins Pen jintegrals defined by equations (11). 
ae a 
i, Mm, Pp, Q) 


positive integers. 
r, 8, k,n} 


A characteristic value phw*b‘/gD. 

Ce) typical element of frequency determinant, defined by 
equations (15). 

ie Kronecker delta defined by equations (15). 

dinm> Pnn parameters defined by equations (17). 

ae parameter defined by equation (20). 


3. The Rayleigh-Ritz method 


In the rectangular plate OA BC, coordinate axes Ox, Oy are taken as 
shown in Fig. 1. The plate is assumed to be rigidly supported against 
transverse displacement around the edges OA, AB, BC, CO, and these 
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edges are elastically restrained against rotation. Opposite edges are taken 
to be restrained to the same degree. If the plate is vibrating harmonically 


ie ee 


Fic. 1. Details of rectangular plate. 


with amplitude w(z,y) and angular frequency w, the maximum strain 
energy is given by 


v= | fa) + Gp +> Se Ge+20—of 5] aar— 
10 f Geashens Sha 


and the maximum kinetic energy is 


a thas ff we aay (2) 


In equation (1) the first term on the right-hand side represents the energy 
stored in the plate, the second term represents the energy transferred 
to the surrounding structure due to the edge rotation. From equations 
(1) and (2) . V 


(ph)29) f f w? dady 
00 





(3) 


The Rayleigh—-Ritz method of solution of equation (3) consists of assuming 
a waveform w(z, y) as a linear series of functions of z and y, and adjusting 
the coefficients in this series so as to minimize equation (3). A set of 
simultaneous equations is thereby obtained, and may be solved for charac- 
teristic values, say A, which determine the critical frequencies of the panel. 
The functions chosen to represent the waveform must satisfy the ‘artificial’ 
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boundary conditions of the problem (5) and, for the panel shown in Fig. 1, 
these boundary conditions may be writtent 


Cw Ep Ow 


ong (y ) u oye 
ow .s 
oy* 


| 
aig 


Along BC (y = 6) 


Along OC (x= 0) w .. aoe 


ex* 
Along AB (x =a) w: 


It will subsequently be shown that a suitable waveform, which satisfies 


the conditions equations (4), may be taken as 
D 
w= §$ $a,,.X0% (5) 
m=1Ln=1 


where X,,, Y,, are the so-called characteristic functions of a beam whose 
ends are elastically restrained against rotation. 


4. The characteristic beam function 


Consider the uniform beam AB, of span 
L, shown in Fig. 2. The ends A and B of 
this beam are rigidly supported and elasti- 
cally restrained against rotation, an equal 
degree of restraint being applied to both 
Fic, 2. Details of uniform beam. ends. The differential equation of vibration 

of the beam may be written 


a‘F, of 








a ~ Ti 
from which 


FL = A. cosh “) +B, sinh “| +C, cos( 7] + sin (7); 


F, is a characteristic beam function and represents the waveform of the 


(7) 


{ Since w = 0 around the rectangular boundary of the panel, then as shown by Leggett 
(6), the expression in the first bracket on the right-hand side of equation (1), which repre- 
sents the strain energy stored in the panel, reduces to 


a b st p 
4D | | (<> + ow) dedy. 
00 


dz* Oy? 


In consequence, numerical values for A may be obtained, which are independent of the 
value of Poisson’s ratio for the panel material. This would not be so, for example, if one of 
the edges of the panel was free. 
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rth mode of vibration of the beam. The constants A,,..., C,, which define 
the modal form, and the coefficient «,, which determines the rth resonant 
frequency, depend on the boundary conditions. For the beam under 
consideration these boundary conditions are 


a GS 
dz? L dx 
@F  é4dF 
é = r == (: ane Gt cd aa 
Meet eS eee 


At z= 0: F, = 0; 
(8) 


Substituting for F, and its derivatives from equation (7) into equation (8), 
then 


A, = —C, = —cot}a,; B, = cot }a,tanh fa, 
a, = —4$£(tan }a,+tanh $a,) 
A, = —C, = tan}a,; B, = —tan }a,coth ja, 
a, = $&(cot 4a,—coth }a,) 


r = l, 3, ete. 
(9) 
r = 2, 4, etc. 
Particular cases of equations (9) arise for a simply supported beam (£ = 0), 
when, after rearrangement, it may be shown that 
A,= B,=C,=0; sna,=0, ic. a, = Pa, 
and for a clamped beam (€ = 00) when 





hin sinh «,—sin «,_ Pees e 


" ~ cosh a,—cos a,’ : ; 


1—cos a, cosh a, = 0. 


Both of these solutions are well known (7). 


Numerical values of «,, A,, and B, up to the sixth mode, computed from 
equations (9) for various values of € are given in Table 1. 


5. Application of the Rayleigh-Ritz method 
If, in equation (5), X,, and Y,, are taken in the form 


X,=A n{cosh (=) — cos(*=)| + By sinh(*=2) + sin(*=*) 


Pi OY) coa(%2)\ 4 B sinh(%Y) 4 sin(*n¥) | 
roo rm 
it is seen that the assumed waveform can be made to satisfy the boundary 
conditions equations (4) term by term and so may legitimately be used 
in the Rayleigh—Ritz analysis. 


5002.45 D 
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It is now convenient to introduce the following notation: 


a b 
| X x ae ki j | Teta dy = M,,., 


0 


a 


dX, , | Pane 


dY,, =e 
“dx dx im? fai J “hed — Kyns 


0 


fy Xn @Y, 
a [x XT fde= En; d f Yi Gh Y= Fw | (11) 
0 


sJliaia” “ae 
0 


l [ d®X, d2X 1[dX, ox) a 





b 
1 f d*Y,, d*Y, ifdY, d?y.7° 
—_ | —-& 2" dy = P en enn 
| dy? dy? © = “*T5] dy |. 
0 





} 


As shown in the Appendix, these integrals have the following numerical 


a* x 


A ae 
Nin = 4(24%— BE +1]+42[B,,+1 ssl 


“mm 


, (12a) 
= 0, ift Am 


= Ey, = Hoy = $08 B3,+1]4+2"[B,—1], if § =m 
ais = oon * (A(By—Day,—Am(B;—l)a;], if i A m; 


i 4 oy 04d. ers | (12b) 


oven’ even 


= 0, ifs +m; m oad, techincal 
odd } 


Similar expressions may be obtained for M,,,, Pi», Kyn, and F,,, by chang- 
ing the subscripts & for i, n for m in equations (12). It is to be noted that 
there is a particular reason why the terms 

1[dX, d*X,,]". dy, @Y,, 

al de dz |.’ b| dy dz? |, 
have been retained in the last of equations (11) ‘aa have not been in- 
corporated directly into N,,, and P,,, respectively. These terms cancel 
with similar terms arising in the expression for the energy stored in the 


structure surrounding the plate and do not, in fact, appear in the final 
solution. 
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Substituting for w(x, y) from equation (5) into equation (3), and using 
the notation of equations (11), we get 


p qa , 
Da ~ —< b 2 
iV ? ° $ 
28 >. 3 Dip Amy b (Nim My» + Lim ken) + *(; Hin, Kyey, - 
~ 44m=1 kn=1 
phabw? << 
t 
—- — S S a;;,4 


im=1kn=1 


mn Lim Mp = 0. (13) 


Differentiating equation (13) with respect to a,,, the resulting equation 
may be expressed in the form 


q , 
> [ ba -A Omn lan _ 0, ( l 4) 
where 
phw*b* 
gD 
i we. see mn 


mm mn 


- 0, if ik ~ mn 


CS 2(2) H K,,, if tk ~ mn 


im 
a 
\ 





b\4 (b\2 
¥( rr 4 4 9 4 fel — 4 
( —_ ° | Xm 1 Xp a M, ni 2) Ts ans if vk = mn 


a 


] 


Equation (14) represents a set of linear simultaneous equations: there will 
be one equation for each of the pq combinations of mn. These equations 
divide into four independent groups as follows: 


(a) those in which 7, m, k, n are all odd, 
(b) those in which 7, m, k, n are all even, 


(c) those in which i, m are odd and k, n are even, 


(d) those in which i, m are even and k, n are odd, 


and each of these groups is associated with one of the four possible com- 
binations of symmetrical and antisymmetrical modes in the coordinate 
directions. 

For a square plate, with all edges restrained to the same degree, groups 
(c) and (d) are essentially the same and lead to identical frequencies. 
Groups (a) and (b), however, can be rearranged to form two subgroups, 
which give independent solutions. The first subgroup includes solutions 
for which @,,,, = @,», Whilst the second subgroup includes solutions for 
which @,,,, = —@nm- The existence of these modes in a square plate 
clamped at all edges is discussed by Warburton (3). 
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6. Solution of equations 


The characteristic values A of the equations (14) may be found from the 
condition that the determinant of the system is zero. However, the labour 
involved in expanding this determinant will generally be prohibitive. By 
using characteristic beam functions for X,, and Y, in equation (5), it is 
found that the diagonal terms in the resulting determinant are much 
greater than the off-diagonal terms. The characteristic values A can there- 
fore be determined by the iterative process developed by Ritz (8) and 
explained in detail by Young (9). 


7. Approximate frequency equation 

As the diagonal terms of the frequency determinant are much greater 
than all other terms, an approximate solution is suggested by representing 
the mnth mode by the mnth term only in the series equation (5). Proceed- 
ing with the energy solution, the mnth resonant frequency may then be 


written anit Fea (Be)[ (Gy ak +a8-+2(5) tam tna? (16) 


where d pe ). sae ree o® | «,,(B?,+-1)+2A,,(B,,—1)] 
mm Lam %m(242,—Bt,+1)+2A4,,(B,,+1) 

K,,, fx, (B2+-1)+2A,(B,—1)] 

M,, %,(2A2—B2+1)+2A,(B,+1) 





(17) 
$nn = 





For any values of the boundary restraints £,, £,, appropriate numerical 
values of «,,,, &»» Pmm: Pnn May be obtained from Table 1. 

It should be noted that, for a simply supported panel, equation (16) 
reduces to the well known exact solution (1) 


vn BQN} 


For a square plate, under equal restraint on all edges, as already stated, 
there exist modes in which a,,,, = a,,, and a,,, = —da,,,, and these occur 
when m is odd/even; n is odd/even, and m #n. The corresponding 
approximate frequency equations may be written 


cinn-nmd = Bau] (5p) Lah + 28+ 2b nO)? 


Wmn+nm) va/ (x) [af +o8+2(Pmm Pant mn) |! 


» <, ieaiiee 
ni } A My», 
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TABLE 2 


Characteristic values vA for square and for rectangular plates 




















































































































































































































| The mode number 
b/a € | I 2 3 4 5 6 
o (| 19°74 49°35 78°96 9870 98°70 
20 | 31°09 64°31 95°85 117°3 116°5 
| 31°16) (64°52) (96°17) (117°8) (116°9) 
I'o 20 35°99 73°41 108-3 131°6 132°3 
| (2611) (73°74) (108-9) (131°7) (132°4) 
| 
° 17°36 41°35 47°47 71°40 S152 I1I4 
20 | 28-21 54°57 61°97 86°85 97°03 123°0 
(28-28) (54°77) (62°17) (87°15) (97°34) (123°3) 
9 co 32°67 62°29 72°76 gd14 109°4 143°5 
(32°78) (62°71) (71°06) (98-66) (109°8) (1441) 
| : | 
-- wee dese td 
° | 16°19 35°14 45°79 64°74 66°72 90°33 
20 25°50 46°02 59°95 79°06 79°24 III*2 
(25°86) (40°17) (60°16) (79°32) (79°59) (11I'5) 
od x 29°08 52°52 65°52 59°40 89°29 124°5 
(29°18) (52°76) (68-80) (89°86) (89°69) (12"°0) 
' ' : 
| ge 
ee 3 oC eas lait bain ae is So 
° | 13°42 24°08 41°55 43°03 53°69 
20 22°30 32°58 50°45 56°97 66°96 
(22°34) (32°68) 50°63) (57°11) (67°17) 
6 oO } 25°g0 37°28 56°93 65°15 75°94 
} (25°97) (37°43) (57°20) (65°39) (76°31) 
} : 
o | 11°45 16°19 24°08 
20 | 20°30 24°15 31°20 
|} (20°33 (24°20) (31°26) 
04 x 22°65 27°81 35°45 
| (23°70) (27°91) (35°56) 
oo 10°26 11°45 13°42 
20 19°38 20°15 21°52 
| (19°39) (20°17) (21°54) 
o"2 co | 22°64 23°45 24°89 
(22°66) (23°49) (24°92) 
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8. Accuracy of solution 


In discussing the accuracy of the present method it is first necessary to 
compare the approximate frequency equation with the more exact series 
solution. For this purpose characteristic values VA have been calculated 
by series solution for plates of side ratios b/a = 1-0, 0-9, 0-8, 0-6, 0-4, and 
0-2, under equal restraint on all four edges (€, = &, = €) and for values 
of € = 20 and . Calculations were based on a 36-term series taking both 
mand n = 1, 2, 3, 4, 5, 6. The results are presented in Table 2. Corre- 
sponding values of vA obtained from the approximate frequency equation 
are given in brackets. For comparison, the exact values of vA for a simply 
supported panel (€ = 0) are also included. It is seen that the approximate 
solution nowhere differs from the series solution by more than 0-7 per cent. 

The series solution overestimates the exact numerical solution to a 
degree which is in general unknown. However, in the case of a square 
plate clamped at all edges Aronszajn (10) has determined lower limits to 
the first ten modes, using a method originated by Weinstein (11). 
Aronszajn actually calculated the value of A/z* but, converting accord- 
ingly, his values for the first six modes of a square clamped plate are 
VA = 35-98, 73°37, 108-2, 131-6, 131-8, 164-4, which differ by only 0-03- 
0-48 per cent from the corresponding upper limits given in Table 2. 
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APPENDIX 


Evaluation of integrals of characteristic functions 


1. The integrals 1, 


4 L 
‘ l "4 1 ‘pee 
z | Fear: i) Prtaes i) rere dx. 
0 0 
L LA L us L 
Fr dx = af FY F, dz = a [Ur Ra | Fi Fi" ar]. 
0 ry 


Now 





. 


0 


But F, = 0 at x = 0 and z = L, whence 


L ms 
| Fae ae [ ere” dx, (21) 
J 


L L 
[ Fr? de = (Fy Fey J P;" F; dz, 
0 
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L L 


and hence | Fy" dx F? dx+[F; F “gd (22) 


| 
It may be shown that 
i. ac* .., 4 
F34— F7"——~ F, F;" = 2(2A}— B?+1). 
Om Xp 
Thus by integration 
L 


| Fac+S | pte tt | Ft Ft” dx = 2L(2A?— B?+1). 
0 = 0 
By substitution from west ions (21) and (22) 
L 
i [ #3 ae }(2A?— B?2+1)- UF; Fy’. (23) 
Stell tr 4ad 
0 


If opposite ends of the beam are restrained to the same degree then, by symmetry 


7 wy 1” pee oat 
(PF, Fs Je=-0 —[F, Fy), heed A,(B, +1). 
Thus equation (23) becomes 
L 
1 = ; 7 A, 
> | Fidx = 24}—B}+1)4+—(B,+)). (24) 
L, a 
0 


2. The integrals 
L 


L - 
L | F.F’ dx; L 


mate 











FP? dz; L | F;"? dz. 
0 6 0 
L L 
We have | RF? de = (FF { FP? de. 
0 6 
But F = 0 atz = Oandz= L 
L ’ 
( FF; dx = — { FP dz. (25) 
” ri) 
L L 
7 pee te L‘ [ r renyL 0rD 
Also, FF, dz = —= | Ur F,'\y— ; aa 
Ny 
0 ny 0 
L L 
009 ” sete at 7 Ld 
Hence Fy” de = (Fy Fr | x, FY?’ de. (26) 
3 “4 
It may be shown that 
L* Made) "9 on pr | 
ae 2F. F; a (Br + 1). 


Thus by integration 

Lf f , 

=| Fv" dx + j F? dx —2 | FF! de = 2% (Bi+1). 
a} L 
0 


By substitution from equations (25) and ae 


L 
, ¥ ” a? ' Pg v7 eit Led 
L | Fade = —L | FF" dx = 3 (B+ Za UF og ol (27) 
0 0 
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If opposite ends of the beam are restrained to the same degree, then, by symmetry 


(F > as hes ° — — (FY ALAA Po Row L = 25 A,( B,— 1). 
Thus equation (27) becomes 
£ I, . 
L | Fede = —1 | RFy de = 2 (B+ 1) +a,4,4B,— 1). 
0 0 
3. The integrals 
L 
1 


0 

Now (a —0G)F, F, = Su, FP,” —F, Fy’ — FF, +F;, F2). 
By integration 

(a$—a pf an dx = [F, F;"—F,F,"—F, Fy +F, Fy}. 

But 


ald é ad 
Fis) = 0; re) = L rs) atx = 0, 
F.,) = 9; Fro = Ep tz=L 
rs) = Ys ne) = — 7 * le atr= tL. 
Therefore 


L 
[ FF, dx = 0. 
D 


L 


L L 
¥ 4 
| FF, dz = & | Fiv F, dx = ral (F." F, ef Fe" F, as]. 
a; 
0 0 


0 
But F, = 0 at x = Oandz = L. Therefore 


L L 
J Ft" Fidx = f F) FY" dx = 0. 
0 0 


L L 
. [| Fe Fide = (F7/ Fy— | Fi Fi dx = 0. 
0 


0 


L 
1 ae ALAA 1 “r , 1 ’ "r 
Thus z |: Fy dx = +(F; Fit = AL 
0 
4. The integrals 
L L 
L | FF; de; L | F, Fy de. 
0 0 
L L 
We have [ BPY de = (F Pag— f FiPFidx (r #8). 
0 0 


But F, = 0 at x = 0 and xz = L; hence 


L L L 
[ RFY dx = | Fr Rae = — | hae 
0 0 


L L 
3 1 F’ Fy” dx 1 we ae 
z | Bhae: LP de; 1 | Peri de (r #8). 
0 


(28) 


(29) 


(30) 


(31) 


(32) 
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L L L 
| 7 pve LI‘ ‘iv pe I er pen eee yere 
F, Fy de = = PY Py dx = | (Fr Fy | Fr" Fy" dz 
0 "9 rt 0 
L 
IA LALA Add Ti alia Diy 
=— |(F; Fy —(F; Fe), + Fy FY dx}, 
Or J 
0 
L 
at 4 ald wee me hope ” 
hence Li FLF, dx =([F, ¥,-—f' fF, i +e , de. 
0 0 
Substituting from equation (32), 
Fs Pi f pep ae — eee Fear Pry 
L| RF? de =L F" F.dx = —L | F{Fidx = ng . (33) 
. ’ J Op — Ag 
0 0 0 


If opposite ends of the beam are restrained to the same degree, then, by symmetry 


?” Ad 
Prez 0 F xa2- L } 





a as r 1, 3, 5, ete. 
Frez 0 o& F adent } 
tpt . =! } yp ; 
tome ne le an 8, 4G ote 
, F ris)z=0 F xe)2~ L } 
Thus 
f FF, dx = 1 frera I f a fee et ee 
, ‘x ae, dx —-L "dz = ris ——; 8 1s —— 
i ts i deri i ee mn even odd’ 
— 40? a’ x5 fou,( B,—1)A,+a,(B,—1)A,] .... odd . odd 
and = Sisenieicatiniammaee ” Ag he ot 6 ip oem, 
of — af even even 
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SUMMARY 


This paper is concerned with the propagation of Love waves in a stratum bounded 
by infinite media with different elastic properties. The symmetrical case is con- 
sidered first and it is shown that there are an infinite number of modes of vibration. 
The properties of the dispersion curves in each mode are derived. These properties 
are illustrated by graphs obtained from numerical computation. It is found that 
the general case is very similar to the case of symmetry, and corresponding results 
are given. The limits between which the phase and group velocities must lie are 
found in each case, and cut-off frequencies and wavelengths are also determined 
for each mode. 


1. Introduction 


SToNELEY (1) has shown that under certain conditions it is possible for 
transverse waves to be transmitted in a stratum of uniform thickness, 
bounded on both sides by infinitely deep materials with different elastic 
properties, without a significant penetration of these waves into the 
infinite materials. Stoneley obtained the secular equations for waves of 
this sort, but his investigation did not yield much further information. 
In this paper it is first assumed that the infinite bounding materials have 
the same elastic properties and with this simplification it is possible to 
obtain results which give us a greater insight into the propagation of these 
waves. In fact, the properties of the phase and group velocities of the waves 
are found to be very similar to transverse waves in a surface layer, which 
were first investigated by Love (2). Once these results have been obtained 
it is not too difficult to extend the investigation to the case where the 
bounding materials possess different elastic properties. 


2. Fundamental equations 


We shall assume that the stratum, consisting of a material labelled 2 
in which the velocity of rotational waves is c,, and the modulus of rigidity 
is #12, occupies the space between the planes z = +h, of a system of 
rectangular cartesian coordinates x, y, z. Materials 1 and 3 extend re- 
spectively from z = h and z = —A to infinity, the corresponding elastic 
constants being ¢,, 14, C3, Hs. If u, v, w are the displacements in the 
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directions x, y, z, then for the propagation of a plane transverse wave in 
the positive x-direction, we may take vu = w = 0, and v in the form 


v, = p,etr- 9 (2.1) 


where 4, is a function of z only. The suffix i = 1, 2, 3 refers to the material 
in which the displacement takes place. The equation to be satisfied by 


the v, is c2V2n, = 2*r,/ét?. 


Thus %,; must satisfy the equation 
d*y,|dz*—(p*—n*/et yy, = 0, 


the solutions of which are circular or hyperbolic functions. We may 
construct the following solution of (2.3), 


y, = Ae", , = Booss,z+Csins,z, and yw, = De**, (2.4) 


where A, B, C, D are constants, and 38,, 8,, 8, are given by 


a = p*—n2/c?, sf = n2/c2—p*, and sf = p?—n?/c3. (2.5) 


Equations (2.4) represent a wave which does not penetrate very deeply 
into the bounding media provided that s, and s, are real and positive. 
Stoneley has shown that if s, is imaginary it is impossible to satisfy the 
boundary conditions. There is no loss of generality in assuming that s, 
is positive and that c, > c,. From (2.5), if ¢ (= n/p) is the phase velocity, 
hen : 
™ (2.6) 
whilst s, < s,. We assume that the tangential displacements and tractions 
are continuous across the planes z = +h. Thus, when z = h, v, = vr, 


> 


ow ov 
and p,— = p,—*, and when z= —h, v,= 15, and ps 
Oz c 53 


~ 


Ov, Ov; 
: = b3- 


2 
>i as Oz Oz 


Substituting the expressions (2.4) in the boundary conditions and eliminat- 
ing A, B, C, D we obtain the secular equation 


m 21 9/¢2 . 27 

8o( 1 8, +73 8,)t - w 2(83 Ta T3 8) 8,)t —8,(7, 8, +73 83) 0, ( f/f ) 

where ¢ = tanhs,, and 7,, 7; are the ratios j1,/y, 3/4. respectively. This 
is essentially the same as the equation obtained by Stoneley, which is 
92 " ¢ 

(83—7, 738, 83)7T' = 8,(7, 8, +7383), 8) 


where 7' = tan 2hs,. 


3. The symmetrical case 
In this section we shall assume that the nedia 1 and 3 possess the same 


elastic properties, so that c,; = c, and yw; = p,. The equation (2.7) now 


reduces to (8,t—7, 8,)(r,8,t-+8,) = 0, (3.1) 
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which means that either 

(i) 7,8, = s,tanhs,, (3.2) 
or (ii) 7,8, = —s, coths,. (3.3) 
The equations (2.7) are now incompatible unless either (i) C = 0 and 
A = D, in which case we arrive at the equation (3.2), or (ii) B = 0 and 
A -D, when we arrive at the equation (3.3). We may therefore 
separate the waves into two kinds, waves (i) which consist of displacements 
which are symmetrical about the plane z = 0, and for which the secular 
equation is (3.2), and waves (ii) for which the secular equation is (3.3) 
and which consist of displacements antisymmetrical about the plane 
z= 0. The equation (3.2) is the same as the secular equation for Love 
waves in a surface layer (2), (3). 

The equations (2.5) now reduce to 


sj = p?—n*/ci, 


8 = n*/c3—p?. (3.4) 


The phase velocity c, p, and n may be obtained from (3.4) in the form 
n? = v,(s8}-+83)c3/(v;—1), 
p? = (v,83+83)/(v,—1), (3.5) 


whilst c? = v,(s}-+83)c3/(v, sf +89), (3.6) 


where v, = c?/c3. By differentiating either (3.2) or (3.3) it can be shown 

that the group velocity V (dn/dp) can be written as 
__ Yal (87-+-89)71 +h (83-+-77 87) leg 
~— [(¥, 8§-+83)7, +48, (83 +75 8p) Je" 

This expression for the group velocity holds for both the symmetrical 
and antisymmetrical cases. 

An examination of equation (3.2) shows that if we take a fixed value 
of s, there is a single solution for s, in each interval of the type 
mn < hs, < (m+4)n, where m is a positive integer, whence corresponding 
values of p, m, c, and V may be found. Thus ¢ and V as functions of p 
depend on the branch of the solution of (3.2) taken, and every branch 
corresponds to a separate mode of vibration. In the same way there 
exists a separate solution of (3.3) in each interval of the type 


(m+4)m < hs, < (m+l1)q. 
Thus each interval }(m—1)m7 < hs, < 4mm corresponds to a separate 
mode of vibration, where m = 1, 2, 3,... and so on. The modes in which 


m is odd consist of symmetric displacements, whilst even m correspond 
to antisymmetric displacements. 


(3.7) 
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Before examining any mode in detail there are a number of general 
conclusions that can be drawn. 
(a2) When hs, tends to the lower limit in each mode, s, tends to zero, 


so that from (3.5) we see that there exist in each mode minimum values 


of p and n given by “page 
i . . p* > 85/(¥,—1), 


and n* > v,(8)c)?/(v;—1), (3.8) 


where hs, = }(m—1)z. Hence the higher modes consist of comparatively 
short waves of high frequency. In particular, if A, the wavelength, is 
greater than 2h/(v,—1)* then this wavelength only occurs in the first two 
modes, whilst only the first (symmetrical) mode can produce a wavelength 
greater than 4h/(v,—1)*. 

(6) Let us consider the behaviour of the phase and group velocities 
when hs, tends to the lower limit }(m—1)z in each mode. From (3.2) or 
(3.3), whichever applies, we can see that s, tends to zero in each case, and 
is very small compared to s,. With the exception of the first mode, which 
is a special case, if we allow s, to become very small in (3.6) then 


¢ = ey[1—(%—1)s¥/s3-+ O(64)], (3.9) 
whilst from (3.7) and (3.9), 


V = c,{1—(v,—1)hs,/7, + O(8})]. (3.10) 
In the first mode, from (3.2), 
8 = hs3+ O(s$), 
whence c = ¢,{1—(v,—1)h?s3/27?+ O(s$)], (3.11) 
and V = e,[1—3(v,—1)h2s3/272+ O(s$)]. (3.12) 

The conclusion to be drawn from (3.9) to (3.12) is that in each mode 
both ¢ and V are equal to c,, when p and n are at their minimum values, 
and that c and V decrease from c, as p increases. 

(c) If hs, tends to the upper limit 4mz in each mode then we see from 
(3.2) or (3.3) that s, tends to infinity and consequently, from (3.5) that 
both p and n tend to infinity, corresponding to very short waves of very 
high frequency. From (3.6), if s, is very large, 

¢ = ef 1+(v,—1)s3/r, 83+ O(s*)], (3.13) 
whilst from (3.7) and (3.13), 

V = cf 1—(v,—1)s3/v, 8? ++ O(s; 4)]. (3.14) 
Remembering our assumption that y, > 1, it is clear from (3.13) and (3.14) 
that as p tends to infinity, c tends to the lower limit c, from above, while 


V tends to c, from below. Since in any one mode both c and V are con- 
tinuous functions of p, V must possess a minimum which is less than c,. 
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(d) From (3.7) we see that V > c3/c, so that on no account can V be 
less than c3/c,. If we fix a value s$ in the first mode, and examine c and V 
in the other modes for s, = s3+-}mz, then t remains unaltered, as will 
the ratio s,/s,. This means that c is independent of m, V decreases steadily 
from its value in the first mode to the limit c3/c. In the limit, if m is very 
large, V = c3/c, except in the immediate neighbourhood of s, = 0. Thus 


c {Mode Mode Mode vw 


en 2 3 
“204 


C2 


























SSE 
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the minimum value of V approaches c3/c, as m becomes very large, and 
the position of this minimum will approach the smallest value of p 
occurring in the mode. 

Fig. | illustrates the dispersion curves of the first 1our modes, with v, 
and +, being taken as 5 and 2 respectively. The graphs, drawn from 
values obtained by numerical computation, illustrate points (a) to (d) 
very well. 

Jeffreys (4) has shown that when there is a stationary value of the 
group velocity then the dispersion is comparatively small for that parti- 
cular value of p. Thus we conclude that though the first signal of trans- 
verse waves will travel with a velocity of c,, and will consist chiefly of 
longer waves, the main part of the vibration will follow, the slowest waves 
travelling with a velocity of c3/c,. In addition, certain frequencies and 
wavelengths, corresponding to the minima of the group velocity, will 
suffer less dispersion, and will predominate at a large distance from the 
disturbance. Further, these predominating waves will travel with v2loci- 
ties between c, and c3/c,, and the important frequencies will be located 
at intervals of roughly c,/4h(v,—1)!. 
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4. The general case 

Here we merely assume that c, > cs > cy, so that s, > s,. In theory 
it is possible to express » in terms of p from the equations (2.5) and (2.7), 
but this is not possible in practice. It is best to eliminate n and p from the 
equations (2.5) to obtain 

v;(v3—1 )si = ¥g(vy—1 )s5 + (vy — ¥3)83, (4.1) 

where v, = c3/c} and v, > v, > 1. We can now proceed by choosing a 
value of s, and determining the corresponding values of 83, s,, p,m, c, and V. 
Before attempting to find particular numerical solutions it is advisable 
to separate the modes of vibration and to try to find some general pro- 
perties of the dispersion curves. The secular equation (2.7) can be con- 
sidered as a quadratic in ¢, with solutions 


(i) t=a+f8, (4.2) 
or (ii) t = a—B, (4.3) 
where x = (7, 73 81 83—8$)/89(7, 81-+73 83), (4.4) 


and B = (s§+-7? 83 5? 


28] +73 83 83 +7] 73 8] 8§)*/(7, 8) +73 85)89. (4.5) 
A glance at (4.4) and (4.5) shows that £ is always greater than |«|, and there- 
fore if hs, lies in an interval m7 < hs, < (m+-4)z, tis positive and the solu- 
tion (i), given by equation (4.2) has to be taken, whilst if (m—4)m < hs, < mz, 
tis negative so that solution (ii), given by equation (4.3), applies. This sug- 
gests that there exists a mode of vibration in each interval of 47, although 
the modes can no longer be separated into symmetrical and antisymmetrical 
ones. We have yet to show that while hs, lies in such an interval there 
exists not more than one solution for s,, 8,, and 8,; this will be done in 
an appendix at the end of this section and for the remainder of this section 
we shall assume that each interval corresponds to a separate mode. 

When s, tends to the upper limit in each mode, 2hs, tends to mz from 
below, so that, from (2.8), 7’ tends to zero through negative values. This 
can only occur if both s, and s, tend to infinity. From (4.1) we see that for 
large 85, 8, = w8,+O0(sz1), (4.6) 
where w is the constant [v3(v;—1)/v,(v,—1)]*. Using this result and the 
equation (2.8) we find that 


T — —Qs, 83+O(s5 %), (4.7) 


where Q = (7, w-+-73)/7,73;. From (2.5) the phase velocity may be written 
in the form 


c? = vg c3(s3+-83)/(vs 8%+-83), (4.8) 
and, if s, is large, 
c = Ce 1+ (vg—1)s3/2v, 82+ O(s5 *)]. (4.9) 
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The group velocity may be obtained in the form 


V = v9(83 83+8)¢3/(vg 83 83+-83)¢, (4.10) 
where s; means ds,/ds,. If we differentiate (4.7) we find that s, = O(s), 
and, using this result, together with (4.9) and (4.10) we find that if s, is 
large, V = cqf1—(vs—1)s3/2r, 83+ O(859)]. (4.11) 
Thus, when s;, and consequently p and nm are large, the phase velocity 
approaches the limit c, from above, whilst the group velocity approaches 
it from below—a result corresponding to the previous section. 

If, as in the previous section, we attempt to make s, tend to the lower 
limit in each mode we find from (2.8) that ifm + 0, both s, and s, tend to 
zero as hs, decreases to 4mz, whilst if m = 0, both s, and 8, are of O(s?) 
as 8, tends to zero. These conclusions are incompatible with the equation 
(4.1), except in case vy, = v3. Excepting this special case, we conclude 
that there is no solution of our equations at the lower end of each interval. 
We now try to find out what happens if we let s, become very small on its 
own. In this case (4.1) can be reduced to 

7,8, = ks,+O(8§), (4.12) 
where k = 7,[(v,—v3)/r,(v3—1)}}, (4.13) 
so that from (2.8), 
, => k+(1+k?*)r, 83/8,-+ O(s?). (4.14) 
We may rewrite (2.8) in the form 
83 = 89(7'8,—T, 8)/73(8g+ Tr, 8,). (4.15) 
Since s, is always positive it follows that if T > 0, 7’ > 7, 8,/s,, and from 
(4.12) it is clear that r,s, > ks,. Hence we see that for all values of 84, 
either 7’ < 0, or T >k. This means that in the mth mode, 
hs, > 4[tan-'k+(m—1)z}. (4.16) 
If, from (2.5), we write n and p in terms of s, and 83, we find that in the 
mth mode the minimum cut-off frequency is given by 
n > vi{tan—k-+-2(m—1)}/2h(vs—1)*, (4.17) 
with a similar expression for p. In particular, the longest possible wave- 
length that can occur is given by 
A > 4hr/|(vs—1)* tan-*k]. (4.18) 
If, in (4.8), we make s, small then we find that 
c = ¢s[1—4(vs—1)s*+- O(s*)], (4.19) 
where s = 8,/8,. By differentiating (4.14) we obtain 
83(1+k) = 2hs,(1+ 77)4+ T—k+ O(s). (4.20) 
E 
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It has been shown that 7' > k in the neighbourhood of s, = 0, so that in 
this neighbourhood s, must be positive and of O(1). From (4.19), (4.20), 
and (4.10) we find that 


V = cg{1—ss3(v,—1)+ O(s?)]. (4.21) 


22, 











0-9 Cot 


Thus, when s, increases from zero, both s, and s,, as well as p and n, will 

increase from their minimum values in each mode, whilst c and V will 

decrease from their maximum value c,, V decreasing much more rapidly 

than c. 

If we differentiate the last two equations in (2.5) with respect to p and 
eliminate V and n, we find that 

de, de, 

V3 83 Je Fo 

dy dp 


Thus v3 838,-+8, > 0, which from (4.10) means that 


= (vs—1)p = 0. 


V > cele > Bley. (4.23) 


From the results in this section it appears that the dispersion curves for 
any one mode look rather like Fig. 2. In fact, the graphs in Fig. 2 have 
been obtained by numerical computation, with v,, vs, 7,, and 7, having the 
assigned values 5, 4, 3, and 2 respectively, in the first mode. In this 
particular case, c; = 2c,, k = 0-77460, tank = 0-66, and p,,,, approxi- 
mately 0-19/h. The numerical work was done assuming values of s,h at 
regular intervals between 0-33 and 7, and then eliminating s, from the 





TRANSVERSE ELASTIC WAVES IN AN INTERNAL STRATUM 51 
equations (4.1) and (2.8). The resulting quartic was solved by Newton’s 
method. Throughout this paper the numerical work has been done by 
means of autocode programmes for the Manchester University digital 
computer. 

In view of the similarity between the results in this section and the 
previous, we can conclude that provided the condition v, > v; > 1 is 
observed, Love waves may be transmitted in a stratum and that at a large 
distance from the disturbance we would expect certain discrete frequencies 
and wavelengths to predominate, corresponding to minima in the group 
velocity. These predominant waves will travel with velocities of between 
c, and c3/c,, although we have shown that the maximum velocity with 
which Love waves can travel is c;. 


APPENDIX 


There remains the question whether we are justified in assuming that there is 
at most one solution of equations (4.1) and (2.8) in every interval 


3(m—1)m < hs, < §mn. 
If we take R = (83—7, 728; 83)/89(7, 8, +7383), (i) 


and we assume that s, is fixed, then 


C dR te 
83(7, 8, +7383)" y (83+7, 73 8, 83)(7, 8, +73 83) — 8,(75 83 + 83)8;, (11) 
2 


where s| = ds,/ds,. A similar expression for dR/ds, in terms of 8; may be obtained 
by assuming s, fixed. From (4.1) we find that s{ < 8,/8, so that, from (ii), if 8, is 
fixed 
cr, > Tz S(T] 8] +83) /83(7, 8, +7383)" > 0. (iii) 
R is a continuous function of s, in the region s, > 0, and when s, is small and 
positive, R is very large and negative. We conclude that the graphs of cot 2hs, 
and of R must intersect once in every interval }(m—1)m < hs, < 4mm, whence 
(4.1) and (2.8) will have at most one solution in such an interval. Similarly, if we 
fix s, we find from (4.1) that 8; < 0, whence dR/ds, > 0, with much the same result 
as above. Finally the question arises as to what happens if we assign a fixed value 
to 8,. From (2.8) we may express 8, as a function of 8, in the form 

7,8, = 8(7'8.—7383)/(8,+ T7383). (iv) 
If we draw the graph of this function, and the function given by (4.1) with s, 
constant, we find that there will be just one intersection for positive s, and 8,, 
provided that either T' < 0, or T > k, whilst if 0 < T < k there will be no such 
intersection. 
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SUMMARY 

An infinite medium of ideal soil contains a single spherical cavity within which 
a slowly increasing pressure is applied. The analysis of stress in the resulting plastic- 
elastic system, in conjunction with the flow rule associated with the Coulomb law 
of failure, leads to an expression for the radial displacement of the soil. It is found, 
in particular, that the radius of the cavity increases indefinitely as the applied 
pressure approaches a finite limiting value. Since the work done at the cavity wall 
also increases without bound this pressure cannot be attained physically. 

When the angle of shearing resistance of the soil is zero the Coulomb law and 
flow rule are, under conditions of radial symmetry, formally identical to the yield 
criteria of Tresca and von Mises and their associated flow rules. A special case of 
the analysis is hence appropriate to the expansion of a spherical cavity in a non- 
hardening metal. Calculations have been carried out for copper and a comparison 
is made with numerical results obtained from an alternative solution due to Hill (1). 
It is shown that the inclusion of work-hardening does not alter the behaviour of the 
medium qualitatively. 

If, on the completion of loading, the pressure exceeds a certain limit further 
plastic flow takes place round the cavity on unloading. Equations are obtained 
from which can be calculated the residual stresses and displacement at any stage 
of the unloading process. 


1. Introduction 


A MEDIUM of ideal soil, unbounded in all directions, contains a single 
spherical cavity. Initially the radius of the cavity is ay and a hydrostatic 
pressure p, acts throughout the soil which is assumed to be homogeneous. 
An additional pressure p is then applied inside the cavity and increased 
sufficiently slowly for dynamical effects to be negligible. This paper is 
concerned with the distribution of stress and displacement in the soil as 
p increases from zero and when unloading subsequently takes place. 

An ideal soil is, by definition, an isotropic elastic—rigid plastic material 
in which the yield condition is obtained from Coulomb’s law of failure. 
When the principal stress components satisfy the inequalities 0; < 0; < o,, 
Shield (2) has shown that the Coulomb law takes the form 


Oy -o, = 2c cos 6 —(o,,-+-0;)sin d, (1) 


where c is the cohesion and ¢ the angle of shearing resistance of the soil. 
(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 1, 1959] 
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1—sin¢’ 1—sin ¢’ 


equation (1) becomes (l+a)o,—o; = Y. (2) 


Putting 


The permutations of the suffixes in equation (2) define six planes in 
principal stress space forming a right hexagonal pyramid. This yield 
surface has also been discussed by Shield. 

When ¢ = 0, then « = 0 and equation (2) reduces to the yield criterion 
of Tresca. This is a special case of some importance as saturated soils 
display no internal friction when deformed in such a way as to preserve 
their natural water content. In the case of radial symmetry the yield 
criteria of Tresca and von Mises and the associated flow rules are formally 
identical. The theory of the expansion of a spherical cavity in Tresca and 

teuss materials therefore emerges as a special case of the analysis which 
follows. 

The properties of an ideal soil are very much simpler than those of 
natural soils and attempts have been made to put the theory of soil 
plasticity on a more realistic basis. Drucker, Gibson, and Henkel (3), for 
instance, have suggested that soils should be treated as work-hardening 
materials. These ideas have not yet been stated in a mathematical form, 
however, and in this paper the resistance of the medium to deformation 
is therefore simply expressed (in equation (1)) in terms of the cohesive 
forces between the constituent particles and the frictional forces set up 
at inter-particle contacts. In the absence of detailed experimental data 
on the behaviour of soils under the conditions envisaged a more elaborate 
description is hardly justified. 

The quasi-static expansion of a spherical cavity in a Reuss material has 
been studied by Hill (1, pp. 103-4). This solution has been applied to the 
determination of bearing capacities by Meyerhof (4), and Skempton, 
Yassin, and Gibson (5) have produced an analysis which incorporates the 
shearing resistance characteristic of ideal soils. The more comprehensive 
treatment presented here is based entirely upon the theory of ideal plastic 
solids and in that respect differs from the earlier work of Skempton et al. 

The initial position of a particle of soil is specified by the spherical polar 
coordinates (r,,@,) measured about the centre of the cavity. Since the 
system has radial symmetry the particle retains the same angular coordi- 
nates throughout the deformation and its current position is given by 
(r,@,%). The principal components of stress and strain everywhere act in 
the directions of r, 6, % increasing. Denoting them by @,, a9, oy; €,, €, €y 


respectively, then G9 = %y, <9 = ¢y- @) 


At the cavity wall, o,(a) = —p—Pp, (4) 
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a being the current radius of the cavity. The other external boundary 


condition is x 
G,(00) = —Ppp. (5) 


The displacement u (= r—r,) is purely radial and the elastic properties 
of the soil are described by Young’s modulus £, and Poisson’s ratio v. 


2. The elastic solution 

As p increases from zero the deformation of the soil is at first entirely 
elastic. Under conditions of radial symmetry the equations connecting 
stress, strain, and displacement are 


du l 
[ 2vaa|, 
ae 
] 
€9 ; A vo,+(1 ~v)a|, (6) 


and the elastic equilibrium equations reduce to the single relation 


, do, ‘7. 


+. 2(0,—99) 0. (7) 


dr 

The solution of equations (6), (7) subject to the boundary conditions (4), 
(5) is ‘a\3 a\3 

o, Po r| ) ; 9 = —pot bel"), (8) 


r 


ui , 2 + v/a’ * (9) 
~ 22 EB ) . . 


From equations (3), (8) it follows that o, < og = oy. Plastic yielding thus 
commences when the state of stress is represented by a point on the edge 
of the yield pyramid given by 


(1+a)og—o, = Y, og = Gy, (10) 


Since (1+ a)og—o, -apy+4(3+-«)p(a/r)> takes its greatest value at 
r = a, yielding starts at the cavity wall when the applied pressure reaches 
the value 2(¥ +app) 


P= — (11) 


3+a 


Both strain components are of greatest magnitude at the cavity wall and 
when p = 7), 








2(1-+-v)(¥ + XP) (l+v)(¥+ XPo) 

eas ’ €9(@) = 7 = . 
(3+a)k (3+a)E 
Soils, like metals, are characterized by the strong inequality Y < £. 
Provided that p, << E the use of the infinitesimal strain theory through- 
out the initial phase of elastic deformation is therefore justified. 


(12) 
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3. The plastic—elastic solution 
> p, a plastic zone is formed around the cavity. Considera- 
tions of symmetry require the outer boundary of this region to be spherical 
and the current value of its radius is denoted by b. 

The stress components in the plastic zone a < r < b satisfy the yield 


condition (10) and the equilibrium equation (7). The solution of these 


When p 


equations is easily shown to be ; 
dl r A 
o, = — + Ar-2all+a), og = —+— po 2al(l+a) (13) 
x a Il+a 


where A is the constant of integration. 
From equations (5) to (7) the stress components in the elastic region 
r > bare 


o, = —p,— Br-, 09 = —Pot $Br-, (14) 
B being a second constant. 

At the junction of the plastic and elastic regions the radial stress com- 
ponent o, must be continuous in order to preserve equilibrium and the 
material on the elastic side of the interface must be about to yield. These 
two internal boundary conditions together imply the continuity of o, and 
og at r = b and determine the constants A, B in the form 





A — 3+ a)(Y + apy) peai+a, Be 2(¥ +p) b3. (15) 
o(3-+ a) 3+0 
The radial stress component in the elastic region already satisfies the 


boundary condition at infinity (5). The other external boundary condi- 
tion, equation (4), gives the relation 


o. yee 

a | 3(1+a)(Y¥+apy) : 

It may readily be verified that 6 > a if and only if p > p,. 
From equations (15), A < 0, B > 0, and it follows from (13), (14) that 

a9 > o, throughout the soil. In the elastic region, moreover, 


(l+a)og—o, <Y for r>b, 





(16) 


the signs of equality holding together. These inequalities establish the 
validity of the solution. 


4. Evaluation of the displacement 
On substituting from equations (14), (15) into (6) the displacement in 
the elastic region, r > b, is found to be 


i oa SEG (17) 


(3+a)E r} 
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The strain components in this region, obtained from equations (6), de- 
crease in magnitude as r increases and attain their greatest absolute values 
at r = b. These maxima are given by the right-hand sides of equations 
(12) so that the strain components at the inner boundary of the elastic 
region remain constant from the initiation of plastic flow in the soil until 
the completion of the loading process. It follows that if py) < # the strains 
in the elastic region are always small enough to warrant the neglect of 
their products and higher powers. 

In the plastic zone indices e, p are used to distinguish the elastic and 
plastic components of the total strains. The plastic strain increments are 
given by the flow rule associated with Coulomb’s law of failure. At all 
points in the plastic zone the state of stress corresponds to a singular point 
on the yield surface. The appropriate form of the plastic potential is 


therefore f = ol +a)og—o,]+(1—y)[(1+a)oy—<,], 


where 0 < y < 1. Differentiation with respect to the stress components 
yields the plastic strain increments 


de? = —d), de} = (1+a)yda, de} = (1l+a)(l1—y)dA, (18) 
A being a non-negative function of r. In general y also depends upon r, 
but the second of the identities (3) requires in the present case that 
e} = «4, and hence y = }. 
Adding to equations (18) the elastic strain increments obtained from 
equations (6), we get 
de, = E-"(do,—2vdag|—d), 
deg = E-| —vdo,+-(1 —v) dog|+-4(1+ a) dA. (19) 
These equations are exact. Eliminating dA, we find 
(1+ a) de,+ 2deg = E-"[(1—2v+-«) do, + 2(1—2v—av) dog]. (20) 
The distributions of stress and strain in the soil at the commencement of 
plastic yielding are obtained from equations (6), (8), (9) by putting p = p,. 
The integral of equation (20) subject to these initial conditions is 
(1+a)e,+2e9 = E-"[(1—2v+a)o,+ 2(1—2v—av)ag+(3+<a)(1 — 2v)po|. 
(21) 
For all a <r <b the plastic strain increments define a normal to the 
yield surface with the direction ratios —2:1+a:1+ a. As these ratios are 
independent of r the effect of further deformation is to translate a given 
normal along the edge of the yield pyramid without rotation. The ex- 


pression of the stress-strain relation of the plastic zone in total form is 
a direct consequence of this property. 
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As r > b, €? and ¢% become vanishingly small and, from continuity con- 
siderations, the elastic strains are also small. In the inner part of the 
plastic zone, however, the total strains become large as p increases and 
the stress-displacement relations (6) cease to be applicable. In order to 
use the stress-strain relation in its total form (21) some definition of finite 
strain must be adopted. The most satisfactory choice is the natural (or 
logarithmic) strain which, in the case of radial symmetry, leads to the 








expressions r™ du 1/du\? 
, ee eee 
2 
a hehe es a (22) 
r, r 2" 
Combining equations (13), (15), (22) with equation (21) there follows 
+0 2a\(1+a 
inf (2) = me p—o(?) eT (23) 
ro dr, r 
where 
ae (o— x)(1— 2v)(Y +-apry) - 3{(3-+ 2a)(1—2v)+-a?}(Y + apy) 
B = exp|- . . , w= . ; 
a(l+a)# a(1+a)(3+a)# 


Since the displacement is continuous at the interface r = 6 the solution 
of equation (23) must satisfy the condition 


1=b = ! ee when r = 6 
(3+a)E 


(see equation (17)). By means of the transformation 


ro (3+a)/(1+a) b\ 2a +a) 
‘ ns (;) : ’ ~ Vr ‘ 





b r 
the solution is found to be 


2a8 y-B+a0/2a (12 reper" (8 ered 
3+a \ (3-+-a)E 





ex{b/r)2a\ +a) 
se f nM +anitagy dy. (24) 


w 


Putting ry = a, r = a, and making use of equation (16), 


2 xB w 8+) 2af 7S. (1 Se en re | 
3+a \ (B+a)E 


a 


aQ 
= J qrareeags dn, (25) 


{(3-+ 2a)(1—2v)-+a2\(¥ apap) 


where Q= Al +atE 
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Equation (25) gives the value of the current cavity radius, a, being 
known. The value of 6 follows from equation (16). The displacement 
r—r, at any radius in the plastic zone can then be determined from 
equation (24). In these calculations the definite integrals must be evaluated 
numerically. 

If the maximum value of p is sufficiently small for the squares and 
higher powers of du/dr and u/r to be negligible throughout the plastic 
zone the displacement can be expressed in the explicit form 





an LLM +o), _ 3((3-+ 2a) —2v)+ oF + “Pel” eae 
rE x(9—a*) E r 

_ 3U—»)(¥ 4 “Pal — 

wan 


When the displacement has been evaluated it is necessary to verify that 
A(r) is non-negative in the range a <r <b. This is the condition that 
an increase in the applied pressure p shall, at no point in the plastic zone, 
result in the performance of a negative amount of plastic work. The 
differential relations (19) can be integrated subject to the conditions pre- 
vailing when p = p,, and since A(r) is then identically zero the integral of 
the second relation is 

L(1-+a)A eg — Ca Pall ey. (26) 
, ak (3-++-a)(1—2v)\r 
The value of ¢g at r = b is given by the second of equations (12) and it 
follows that A(b) = 0 as required. Differentiating equation (26) with 
respect to r, 
Mta)r dn deo 6(1—2v—av)(¥ ee - 
? dr dr (l+a)(3+a)Ez r 





If a < (1—2yv)/v then dA/dr < 0 and A is non-negative in the plastic zone. 
If « > (1—2yv)/v it must be verified by direct calculation from equation 
(26) that A(r) > 0 fora <r <b. 


5. Application of the theory of loading to metals 

In the limit «0 the equations derived in the preceding sections 
describe the response of a Tresca (or Reuss) material, initially at uniform 
hydrostatic pressure py, to internal static loading at the face of a spherical 
cavity. The treatment of this special case is simplified by defining the 


—— 8 = 2(1—2)¥/E, «= }(1+»)¥/E. (27) 
For most metals these numbers are of order 10-3, 


Plastic flow commences when the applied pressure reaches the value 3Y, 
where Y is the yield stress of the material in uniaxial compression or 
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tension. When p exceeds this value the components of stress in the plastic 
zone a <r <b are 


o, = —pyo—RY 4+ 2Y In(r/b), og = —pPot+4hY+2Y In(r/d), 
and the radii a, 6 satisfy the relation 
b 2 : 
= exp( Jp -3): (28) 


Using the method of section 4 to evaluate the displacement in the plastic 
zone the differential equation corresponding to (23) is found to be 


r\? dr re r - 
in| (7) i = 381m. (29) 


The continuity requirement on the displacement at the interface r = 6 is, 
from equation (17), 
ry = by = (l—e)b_ ~whenr = b, 


and the solution of equation (29) subject to this condition is 


(;) peter ~8){(j)' -a—e"]. (30) 


At the cavity wall, r = a when r, = a, and, making use of equation (28), 


exp| l —a)(1 _ oy) | —1] = (1 ~8)|(%2) exp(1 _ | —(1—er|. (31) 


The quotient a/ay is easily determined from this equation as a function of 
p. The value of } is then given by equation (28) and r, may finally be 
calculated from equation (30) when r is specified. 

It will be seen from (31) that a/a, is infinite at a finite pressure p, 


gran by Po = §Y[1—(1—8) nf{1 —(1—8)(1 69}. (32) 
This behaviour is properly understood in terms of the integral 





Pp 
W = $n | p da’, 
ra 


which represents the total work done in raising the pressure on the cavity 
wall from p, to p. Substituting for a* from equation (31) it is not difficult 
to prove that W—+oo as p-p,. The limiting pressure p,, is hence 
physically unattainable. 

Values of p,,/Y for various combinations of Y/Z and v are contained in 
Table 1. These results indicate that neither metals nor soils can, in equi- 
librium, support internal pressures of more than a few times the com- 
pressive yield stress. They also show that p, varies only slightly as the 
value of Poisson’s ratio changes. All the other calculations presented in 
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this paper refer to the expansion of a cavity in a medium of heavily cold- 
worked copper. It is assumed that no further work-hardening occurs. 
The appropriate constants, Y/H = 2-19x 10-4, v = 0-34, are taken from 
a paper by Bishop, Hill, and Mott (6). The variation of a/a, with p/Y for 


TABLE | 
The variation of p,.}¥ with Y/E and v (a = 0) 





o’3 





37205 


4°78! 
6°313 











7 Ost 





copper is given in the first two columns of Table 2 and represented 
graphically in Fig. 1. It is seen that a/ay increases extremely sharply 
as p approaches the limiting value p,. On repeating these calculations 
under the assumption of incompressibility (i.e. vy = 4) no substantial 
differences were found. The values of a/a, are somewhat smaller, as would 
be expected, and the limit p,./Y is raised to 4-480, 


TABLE 2 


The variation of aja, and b/a with p|Y for cold-worked copper 





ala 





Present Hill's 


solution solution 





"0010 ‘0010 
*OO17 0017 
+0040 ‘0040 
"00904 "0095 
‘0217 ‘0215 
"0504 oso6 
"1242 "1250 
"3983 "4036 





4°301 4°296 














As mentioned in section 1, an alternative treatment of the expansion 
of a spherical cavity in a metallic medium has been given by Hill (1, pp. 
103-4). An earlier analysis by Bishop, Hill, and Mott (6) neglects the 
compressibility of the material, but is then able to include the effects of 
work-hardening in an approximate way. The availability of these related 
investigations prompts the interesting question of the exten! to which 
the behaviour of the medium depends upon the particular mathematical 
methods and physical approximations employed in the analysis. In the 
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remainder of this section the results of Hill et al. which bear upon this 
question are briefly recapitulated and somewhat extended. 


2-07 


1-9 





8 
he] 





ig men nn nnn nnn nnn nnn nnn nnn 


10 ; : 


00 05 10 15 20 25 30 355 40 + 
—/p/Y}—> 
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The mathematical method of Hill differs from that used here in its 
treatment of the displacement in the plastic zone. Retaining the incre- 
mental form of the stress-strain relation and then evaluating each increment 
following the motion, by using the interface radius b as a measure of the 
progress of the deformation, Hill arrives at the differential equation (in 
the present notation) 


da b\? fa 
5 = (6+39(7) -3(5). (33) 


At the commencement of plastic-elastic deformation, 


a= 6b =a,(l—e)", 
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(from equations (9), (11)). The solution of equation (33) subject to this 
condition is 


: (l—e)“ex (; p 1 —3e 1a+5)~ 
i iiiate. 3 2Y)| (1+8)exp{1—(3p/2¥)}— (5+ 3e) 


(34) 
and it follows that a/a, is infinite at the finite pressure 


Pa [1 + in( 5) (35) 


Results calculated for cold-worked copper from equations (34), (35) are 
given in the final column of Table 2. They are in good agreement with the 
values obtained from equations (31), (32). 

For the special case « = 0, appropriate to metals, the most serious 
physical simplification made in the present analysis is the neglect of 
strain-hardening. The inclusion of this effect is achieved by writing the 
yield condition in the form 


og—0, = Y+H( ( de), og = Fy, 


where de? = [§(deP)*+ §(deG)?]* 


is the increment of equivalent plastic strain. The stress analysis proceeds 
exactly as before and provides the result 


b 
p= aY(1 +3In “) + 2 | H( [ den)" (36) 
a J P r 


which reduces to equation (28) when the work-hardening function H is 
identically zero. In order to find a an additional relation between 6 and 


a is required, This is easily found if the medium is assumed to be incom- 
pressible. The condition of mass conservation then takes the form 


—_r — ai—ai = 6'—b? = 763, (37) 


Y\3 
where 1=1 (1-3 = 1—(l—e)*, 


the last equality following from equation (17) with a= 0, v=}. As 
stated above, the neglect of compressibility does not by itself produce 
wide deviations from the exact solution. By means of the result 


deP == 2 €9!\; 


which may readily be shown to follow from the assumption that the plastic 
strain components are much larger than the corresponding elastic strains, 
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equations (36), (37) can be combined in the form 


p = $Y(1+Inp—Inz)+4 f afamo(,4))%. 


a,\3 
where == }—{-9} , 
cm 


-1 
It follows that i os | ¥+H(2m*)] > 0, 


whence a/a, increases monotonically with p. Thus if p is finite when 
#. = 1 the behaviour illustrated in Fig. 1 is reproduced in a work-hardening 
medium. Putting » = 1 in equation (38) and replacing p by p,, an 
obvious transformation of the integral leads to the equation 


H(n) dy 


= 3Y(1—Inr Bet.) teak: SR 
ttt iTine ihe | exp(}m)—1 


—2in(1—e) 


(39) 


The integral on the right clearly converges for a wide class of functions 
H(n). In actual fact the work-hardening properties of most metals can be 
accurately represented by a polynomial expression of the type 


nl fa) = $m fa 


n=1 


The corresponding form of equation (39) is 


N 
! ad 
» = #Y¥(1—Inr)+ H,, i) POS. icine SOS 
oa eee 2 exp) —1 
—2In(1—«) 
and since r< 1, the lower limit of each integral can, with sufficient 
accuracy, be replaced by zero. It then follows (see (7), pp. 265-6) that 


N 
Po = §Y¥(1—In7r)+ 2, (%)"+4 n! C(n+-1), (40) 
¢ being the Riemann zeta function. In the case of linear work-hardening, 
N = 1 and since {(2) = }n?, 
Po = §Y(1—In7)+-30°A,, 


a result due to Bishop, Hill, and Mott (6). From these results it can 
safely be concluded that non-hardening and work-hardening media re- 
spond in qualitatively the same way to pressure applied in an internal 
cavity. 
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6. Unloading 

Returning to the case of an ideal soil and assuming that the pressure 
in the cavity has been raised monotonically to the value p (< p,.), con- 
sider now the effect on the components of stress and displacement of a 
monotonic unloading to the pressure (l—zx)p (0 <x <1). The con- 
figuration of the system on the completion of loading is used as a reference 
state from which the variables of the unloading process are measured. 
Deviations from this state are denoted by a single prime and absolute 
quantities in a partially unloaded state by a double prime. Thus 


” , 


' , ” , ‘ ” ' 
Oo, = 0, + 9,, Og = 99+ %, u=—Ust. (41) 


Consistency with the notation used in previous sections is achieved by 
using unprimed symbols to describe the reference state. 

The work expended on plastic deformation during the enlargement of 
the cavity is not recoverable on unloading. Thus when the pressure is 
reduced from its peak value p the medium unloads elastically and 


; a’\3 ‘ 1 a”\8 
sc sonal (4) ’ 09 = —4er(5,) ’ 


1+v/a’\3 
uo =r—r” = rp —_i-—] rv’, 42 
oF E (" (42) 


” 


these quantities being the elastic components of stress and displacement 
produced by the application of an additional pressure —p at the cavity 
wall r”° = a". Had no plastic deformation occurred during loading (i.e. 
if p < p,), the residual stresses and displacement would be obtained by 
adding these expressions to equations (8), (9). When completely unloaded 
(i.e. when x = 1) the system would then have reverted to its initial state. 
If p>p,, however, the residual stresses are, in the plastic zone 


sar =e. 
of —Y_ 31 +ay(¥4 OO a(t) 





O x(3-+-a) r rv")? 
‘4 « ri 2a/(1+ a) "\3 

wey ee ee —t2p(5) é (43) 
x a(3-+-a) \r r 


and, in the elastic region r” > b”, 


2(¥+ap,)(b\® , _ (a”\? 
ae Pol +2n(%;) ’ 


3+-a 
7 Y+ b 3 3 
o9 = —Pot P| —Ie(5) . (44) 


As x increases o; increases and o” decreases. Thus (1+a)og—o; < Y 
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throughout the soil, but if p is sufficiently large of may eventually exceed 
og. It is then necessary to examine the possibility that the yield condition 

(l+a)o;—og = Y (45) 
is satisfied at some level in the soil before the completion of unloading. 
Expressing equations (43) in terms of r” and a” by means of equations 
(16) and (42), there follows 


(1+-a)o> e: 
(3+a) (14a) 
= r+(5 ‘) | r0+20)— F207 +ap+ape)( a x 


[ort yiee PET 


for a” <r" <b". As the expression on the right is a decreasing function 
of r", vielding recommences at the cavity wall when x reaches the value 


- 2(2+-a)(Y +-ap+ apy) 
. (1+-a)(3+ 2a)p 
It can easily be verified, by means of equations (44) together with (16) 
and (42), that the yield criterion (45) cannot be satisfied at any level in 
the elastic region when x < 2,. If p > p,., where 


2(2+a)(¥ +capp) 

3+a ‘ 
then 0 < 2, < 1, and a secondary plastic zone must be formed round the 
cavity during the course of unloading. If p, < p < pg, then the medium 
can be completely unloaded without the occurrence of further plastic flow. 
For values of p within this range the residual stresses are given by equa- 
tions (43), (44) for 0 <a <1, and the residual displacement r”—r, is 
found by combining equations (24), (25), and (17) with (42). 

Consider now the situation when p > p, and x > z,. Let c” be the 
outer radius of the secondary plastic zone and assume, for the present, 
that ce” < 6”. In the range a” <r” <c” the residual stress components 
o”, og satisfy the equilibrium equation (7) and the yield criterion (45). 
The general solution of these equations is 





(46) 





P, = (47) 


of = Ate, ope Ati tayrt, (48) 
a a 


Since the remainder of the soil has unloaded elastically, 
o, = B’r’s, =o = — 3 B’r’ (49) 
for r” >". The conditions of equilibrium and incipient yielding at r” = c” 


imply that of and o% are continuous across this interface. These conditions 
5092.45 F 
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determine the two disposable constants A”, B” in the form 





A’ ei aes ee ore, B’ = Pa 
a(3-+ 2a) \c (3+a)(3+2a) \c 
(50) 


The condition at the cavity wall, 
ofa") = —py—(1—2)p, (51) 

provides the additional relation 

a”\24/q\20A+e) — (1+a)(3+2a){¥ +-ap(1—x)+ apo} 

(F) () Pe (3+-a)(Y+ap+apo) 

The unloading displacement in the region r” > c” is 


Cc 
ut — —B'ty) _ _ 32 +a)(1+v)(¥+ap9) (62+ /(c"\*_. (59) 
~ Qkr® — (84+a)(34+2a)F \c ’ START, 








” 


= 
and, in particular, 
= : 3(2-+a)(1-+-v)(¥ +-apy) (b\ 200+ 
e=c"—u'(c”) = cj} 1+ 3 — rv\(d tee ; 
(B+a)(3+4+2a)E \c 
The right-hand side of equation (52) is less than unity if and only if 
x > x, and since x is known to satisfy this inequality, 


(3. 2a a 2a/(1+a) , 
c Cc 


As a" <c” if and only if a <c, it follows that a” <c”. Due to the 
occurrence of secondary flow at the cavity wall the value of u’(a”) = a—a" 
exceeds the elastic unloading displacement at r” = a” which would be 
given by equation (53). Hence a/c > a”/c” and, using equations (16), (52), 
( 2a b 2a/(1+a) i a” 2a b 2a/(1+a) = (3+ 2a){¥+ap(1 —#)+-apo} > 1 
c} \e “\e") \e 3(¥ +ap,) a 
Hence b > c and a” <c” < 6". These inequalities confirm that a self- 


consistent solution for the stresses follows from the assumption that c” < 6” 
made earlier. 





(54) 





There remains the evaluation of the displacement in the secondary 
plastic zone. The plastic potential associated with the yield condition (45), 


f" = y'[(1+a)o,—9]+(1—y")[(1+-a)o,— oy], 
generates the plastic strain increments 
de’? = (1l+a)d)’, dey? = —y" dr’, dei,” = —(l—y’)dd’. 


Applying the symmetry requirement that y” = 4 and adding the elastic 
strain increments, 


de? = E-"(dol—2vdog|+(1+a)dd’, 
deg = E- —vdol+(1—v)dog|—} dd’. (55) 
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Eliminating d)’, 
de*+2(1-+0)deg = E-{{1—2v(1+a)}do%+2{(1—2v)+-a(1 —v)} dos]. 
(56) 
The differential relations (55), (56) are exact and can be integrated subject 
to the conditions attending the commencement of secondary plastic flow. 


The latter are obtained from equations (42), (43) by setting 2 = x, and 
the total stress-strain relation is found to be 


+a)eg = | 0 —24(1-+a)jof +2{(1— 20) +a(1—v)}05— 


(3+2a)(1—2v)¥ | 3{3(1-+a)(1—2v)—2a%}(¥ a (57) 
- 





x a(3+ a) 


The logarithmic strain components measured from the fully loaded state 


are r” 
o= in(“). (58) 

r 
Entering these expressions into equation (57) and substituting for 07, o% 


from equations (48), (50) there follows the differential equation 


int r”\ 2+ "dr" a __ 3{(3+ 4a)(1 — 2v)4-2a7(1 —v)}(¥ + XPo) b\ 22 +a) /p”\ 20 
5 * "S Geoteeaity a(3-+2a)H c c” 





c 





a(3-+a)E 
This is to be integrated under the boundary condition, r” = c” when 
r =c; ¢, ce" being connected by equation (54). By means of the trans- 


formation c\ 2a(l+a) r”\20 
=i y= a(5) 


_ 3{3(1+-a)(1—2v)— 2a*v}(Y + apy) b 2a\(1+a) 
a a(3+-a)E c , 


4 A at Ser oe) oye (59) 


where 





xX = 





e 3{(3+4 ee ae 


a(3+2a)E 


the solution is obtained in the form 
x(e/r)tatr* a) 


La) £-GB+aX1+2a) 2apt dé +- 


«x 


c 


x(rie) 


a" 3+2a . pste¥ ‘ 
4. (-) K~(1+0KX3+ aN2a3+8a)80 | E® 2a dE = 0. (60) 


x 
Putting r = a, r" = a” and making use of equations (16), (52), and (54) 
an equation for 6/c is obtained. Since the value of 6 is known from the 
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loading solution this equation specifies c, and c” follows from equation (54) 
so that r” can then be determined as a function of r from equation (60). 
After evaluating the displacement it must be verified that A” is non- 
negative in the range a” <r” <c”. The second of the differential 
relations (55) has the integral 
1)” ies 3() + apo) - 
5 6 EB AN 
2 


9,, 9 2a(1+a) 9,,) 1 —" 2a/(1+a)/p”\ 2a 
- | 1—2: (°) (1—2y) e—() (=) | (61) 
3+ a r 3+ 2a Cc Cc 


Since the tangential strain component at r” = c” is u’(c”)/c” it follows from 
equation (53) that A"(c”) = 0. Differentiating equation (61) with respect 
to r”, we obtain 

1M dr" be de @ 6(Y 4 XPo) 

= dr dr E 


s/E\ 4 


1—2v—av [b\2?*%4+%r” dr (1—2yv)- oe p"\ 20 
: (1-+-a)(3-4 (5) r dr" ’ 3+ 2a ) c” : 


Since —e% is a decreasing function of r”, the condition « < (1—2y)/y is 
sufficient for A” to be non-negative in the secondary plastic zone. If 
x > (1—2v)/v the necessary verification must be carried out numerically 


by means of equation (61). 


7. Application of the theory of unloading to metals 

For the case of a metal (i.e. when a = 0), plastic flow occurs during 
unloading when x > 4Y/3p. The stress components in the secondary 
plastic zone a” <r” < 6” are then 


; br” >, [or” 
o. »-+- 2) —ZY in : op = —Py—hY —2Y In|—}, 
r Pots ce” 0 Po 5 cc 
where a”, c” satisfy the relation 
aa” 2 2p . 
= = exp|—-— P , (62) 
ce 3 2Y 


The method used in section 6 to evaluate the displacement in the secondary 
plastic zone leads in this case to the differential equation 


in| ( 5 | > —3d in(), (63) 
rj} dr- ce 


use being made in this section of the definitions (27). 
The solution which satisfies the condition r” = c” when r = ¢, is 


a—a(2)[(G)""- | a a+ay{(")" es i), (64) 
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where, from equation (54), 


Cn (1 (65) 


, 


c 


” 


Putting r” = a",r = a in (64) and making use of equations (62), (65), there 


follows 
Fe 3rp Cc 3(1+8) 1+ 5 a 31-8) 
xp] (1-+8)(2— 2"? }1/° —1 = (142) — <I (- —1}. (66 
alos a7)l) —'Oter St) — I] 
Having determined c from this result, a being known from the loading 
solution, the values of c” and a” follow from equations (62) and (65). 
r” can then be evaluated from equation (64) when r is specified. The 


a, ” 


displacement in the region r” > c” is obtained from equation (53) with 


TABLE 3 
The variation of a"/a and c"/a" with p/Y for cold-worked copper after 
complete unloading 





ply a’ la 
4 09980 
I°5 09978 
2°0 0°9966 
2°5 09948 
30 "9920 
3°5 09879 
40 09817 














TABLE 4 


The variation of r/rg with r/a in cold-worked copper after loading to p/Y = 4 
and the corresponding values of r"/r, and r"/a” after complete unloading 





After loading After complete unloading 





Description of Description of 
material r/ro r’ |r material 





Cavity wall r =a ° 1-3983 ‘ 1-3727 Cavity wallr” = a” 

Plastic “2! 1-1381 +22 1-1333 Secondary plastic 
a 1-0700 , 1-0684 a 
1-0411 . 1-0404 a 


1-0275 Interface r” = ¢ 
1-0263 2- 1-0259 Plastic 
1-0068 , 1-0066 sl 
1-0025 . 1-0024 o 
1-0012 - 1-0011 ie 
Interface r = b 1-00098 : 1-00088 | Interface r” = 6” 
Elastic 5 1-00067 . 1-00064 | Elastic 
9° < 1-00042 ° 1-00040 
1-00028 S 1-00027 
1-00020 . 1-00019 
1-00015 . 1-00014 
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«a = 0. When” is known as a function of r throughout the medium appeal 
to the loading solution gives the residual displacement in terms of ry. 
Some further calculations for a medium of heavily cold-worked copper 
with the properties listed in section 5 are presented in Tables 3 and 4. 
Table 3 gives values of a”/a and c/a” after complete unloading from a 


number of peak pressures in the range {4¥ <p <4Y. These results 


1-457 
1-40 -------------------- 


1°35 


T 





Plastic zone Elastic region 


1:°05-+ 








1:00 ——— er 
' 


demonstrate the very small extent to which the cavity radius is reduced 
when the pressure is returned to its initial value. Table 4 contains values 
of r/r, at distances between 1 and 10 cavity radii after loading to a peak 
pressure of 4Y together with the corresponding values of r’/r, after 
complete unloading. The distribution of r/r, on the completion of load- 
ing is shown graphically in Fig. 2. The displacement is seen to decrease 
extremely rapidly near the cavity wall and this implies that finite strains 
are produced only within the innermost half of the plastic zone. 
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SUMMARY 

This paper applies the theory of sectionally holomorphic functions, developed by 
N. I. Muskhelishvili and others, to some mixed boundary-value problems of aeolo- 
tropic thin plate theory. The mixed boundary conditions are derived from plates 
having boundaries which are partly clamped, and partly free or subjected to specified 
bending moment and shear. Each problem is reduced to the solution of a series of 
equations relating the boundary values of sectionally holomorphic functions on 
opposite sides of one or more ares; these are problems solved in (2). Details of the 
solution are given for the problem of a plate in the form of the upper half-plane, 
subjected to constant bending moment and shear on a segment of the real axis, and 
clamped along the remainder. 


1. Introduction 

IN a previous paper (1) the author considered some mixed boundary-value 
problems in isotropic thin plate theory; in this paper some of the corre- 
sponding problems are studied for aeolotropic plates. It is shown that 
the mixed boundary-value problems may be reduced to the solution of 
two non-homogeneous Hilbert problems ((2), part IV) whose general 
solutions are known, and that in this case the solutions of special problems 
may be also simplified by using contour integration. In the greater part 
of this paper the plates are considered to have the most general type of 
digonal symmetry (i.e. elastic symmetry about the midplane of the plate) 
but even so the complex variable analysis leads to a solution reasonable 
free from burdensome algebra. 

The plates considered occupy one of the following regions: 

(i) the upper half plane; 

(ii) the infinite region bounded internally by a circle; 

(iii) the infinite region bounded internally by an ellipse. 

(Note that the boundary-value problem for an aeolotropic circular disk 
is similar to the problem of an isotropic elliptic disk, and cannot be dealt 
with by the present method.) The plates are clamped along a part L, of 
their boundary L and are free, or have specified bending moment and shear, 
on L,, the remainder of L. There are two main types of mixed boundary- 
value problems, namely, 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 1, 1959) 
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(a) L, free and the plate subjected to an isolated normal load at an 
interior point; 
(6) L, subjected to specified bending moment and shear. 


2. Fundamental equations 
Consider first the case when the material has the most general type of 
digonal elastic symmetry, the stress-strain relations containing thirteen 
elastic constants ((3), 158). Referred to rectangular cartesian axes 
O(x,y, Z) with the origin in the mid-plane of the plate (bounded by 
Z = +h), the transverse displacement w of the mid-plane satisfies the 
fourth order differential re 
dw 3p 
saont Sia a oz! C2 
where a, and a, are real constants, a, and a, are complex constants, 


z = x+1y, and pis obtained from the boundary conditions over the faces 
of the plate, viz. 


wor ‘42a 2a . sgt (2m +a,)— = +435 (2.1) 


= 


tz = 0 = 9 = 22 over Z = —h, rz = 0 = yz, 22 = —pover Z=h. 
(2.2) 
Equation (2.1) is given in a different notation by Green and Zerna (3): 
it may be noted that when the plate is orthotropic the constants a, and 
a, are real. If the equation 
y+ 2a, y+ (2a,+a4)y*+ 24,7 +4; y* = 0 (2.3) 
has distinct roots y,, y, satisfying |y,| <1, |y,| < 1, then the general 
solution of (2.1) is 
w= WytD,(2,)+2,(2,)+0,(z,)+Q,(z2), (2.4) 
where z; = z+y,;2 (j = 1,2), Q,(z;) is an analytic function of z; in the 
region of the z,;-plane which corresponds to the region of the z-plane 
occupied by elastic material, and w, is a particular solution of (2.1). If 
(2.3) has a repeated root y satisfying 0 < |y| < 1, then the general solution 
OF 21) Iw = 9+ 5, le) +2, 2G) +0(4) +0) (2.5) 
where z, = z+y2 and Q(z,), w(z,) are analytic functions of z,. It may be 
shown (Gladwell (4)) that, in all the problems treated for the general type 
of aeolotropy, the solution of the problem for this special type may be 


obtained by a simple modification of the method used in the isotropic 
case. 


3. Boundary conditions 
The boundary conditions along L,, the clamped part of the boundary, 


_ w= 0 = dw/én on L,, (3.1) 













74 G. M. L. GLADWELL 


which may be written 


w= 0= dw/éz on Jy. (3.2) 
Using the notation of Stevenson (5), the boundary conditions on L, are 
ns = M, nz—énn|/és = S_ on Lz, (3.3) 


where M = 0 = S when J, is free. 

Equation (2.4) allows the boundary conditions (3.2) to be written in 
terms of the complex potentials; it is desired to transform (3.3) in the 
same way. If F and K are respestively the resultant shearing force along 
an arc of Z and the complex moment about the origin of the stresses along 
that arc, both measured as per unit thickness of the plate, then it is shown 
in (1) that (3.3) is equivalent to 


re (K+icF) = its= Mon Ly, (3.4) 
on Pe yee ~ Onn ‘ 
im— —(K+izF) = i7z—-—-=S on I,. (3.5) 
és dz és 


These equations show that specification of bending moment and shear 
along an are enables the determination there of K+izF apart from an 
arbitrary combination of the form iCz+-e, where C is a real constant and 
« is an imaginary constant. In particular the boundary condition for a 
free boundary may be written 


K+izF =iCz+e on Z,. (3.6) 
It may easily be shown ((4) ef. (3), 252) that 
y . 8h? , ee | , . ro ee ~ 
K+izF = - 5 {Aq 4 (24) pty 24 (24) + Ag Q9(Zp) + He 229(22)}, (3.7) 
where 
Ay = (a, yf +4, y;+43)/7;, By = A4+4,7;+0377 (j = 1,2). (3.8) 


By using (2.4) and (3.7), the boundary conditions (3.1) and (3.3) may now 
be written 


Q,(2,) +04 (2) +-2,(z)+0,(z,) = 0 on Ly, (3.9) 
OY (2) +7 Q(z) +04 (22) + 7.224(z_) = 0 on Ly, (3.10) 
Ay 4 (24) + poy Q4 (1) FAQ DG (Z2) + pg. Qh (Z_) = M(z,Z)+iCz+e on Ly, 
(3.11) 
where M(z, 2) satisfies the equations 
Sh? | 

T(guea) =M on J,, (3.12) 

8h?. (ad 
im(s 7, Mte3) =S on J,. (3.13) 

: és dz 


In particular, if L, is free (case (a)) and QF(z;) is the complex potential 
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for the corresponding problem when the plate is clamped along the whole 
of L, then writing 


2,25) = QF(z;)+OF(z;) (7 = 1,2), (3.14) 
the boundary conditions become 
OY (24) + F1 QY (21) +08’ (22) +-72.02'(z2) = 9 on Ly, (3.15) 
Ay 29 (z,) + py eq) + Ag 28’ (2) +149 eq) = M(z,Z)+iCz+e on L,, 


(3.16) 
where o - AT. e - 587 
OF (2) +7, QF (21) + OF (22) +72 3'(z2) = 9 on L, (3.17) 
A, OF (24) py QP" (24) Ag Q3" (2) + tg. Q3" (2) = —M(z,2) on Ly. 
(3.18) 


As in the isotropic case, it is assumed here also that the displacement w 
is at most O(|z|) at infinity, and inspection of (2.4) shows that this implies 
Gat Q,(2,) = iHj2?+O(2,), (3.19) 
where the H; satisfy 

H,—H,7i+H,— Hy = 0= Ay,—H,y,+Ayy2—Hey2- (3.20) 
Since in case (a) it may be shown that the (known) Q}F(z,;) satisfy such 


a condition, it follows that the Q?(z,) also satisfy the condition. 
Further, it is easy to show that if 


A,—A,7,+A,—Asy2 = 0 = B,—B,+ B,—B,, (3.21) 
then the complex potentials 
2 
give w= 0, K+izF = om (iz+«) (3.23) 


for some real C and complex «. For these reasons the terms iCz+e in 
(3.11) and (3.16) may be neglected, and the functions ’(z,;) in case (a), 
and '(z;) in case (b), may be taken to be holomorphic functions of z,; in 
the region occupied by the material, including the point at infinity. 


4. Mapping of the regions 

Denote by 7'+ the region occupied by the material in the z-plane, and 
let T’- be the complement of 7'*++ L. The mapping functions z; = 2+-y,Z 
(j = 1,2) transform 7'+, 7, and L in the z-plane into 7}, T;7, and A, 
respectively, in the z,-planes. It is desirable to consider boundary condi- 
tions along these new boundaries A;. In order to effect this, the regions 
Tj are mapped on to either the upper half planes (in case (i)) or the 
exteriors of circles (in cases (ii) and (iii)) in ¢,-planes. Let t; represent an 
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admissible parameter on the boundary in the ¢,-plane; identical values ¢ 
of ¢, and t, can be made to correspond to the same value of the parameter 
z on L by a proper choice of the mapping functions z,(Z;). The ¢,-planes 
may thus be considered as superimposed, and the regions corresponding 
to T7, T;, and A; may be denoted by S+, S-, and C respectively: C, and 
C, are used to denote the parts of the boundary C which correspond to 
L, and L,. 

The mapping functions are now given for the three cases: 

(i) Z is the real axis, and the mapping function is ((3), 332) 


z,(C;) = 2+ V;2 _ (1+-y,)¢;. (4.1) 
It is seen that if z = x+y then 
“= r+(i iy (4.2) 
l +Y¥; 


and therefore, on the real axis {, = (, = 2. 
(ii) If LZ is the circle |z| = a then ((3), 348) 
z;(C;) i z+y;2 —_ C;+4 y;@*/C;. (4.3) 
Since | y;| < 1, the transformation has no singularities (i.e. points at which 
2i(,) = 0) in |f,| >a. 
(iii) If L is the ellipse which is given parametrically by ((3), 364) 


2,(f;) = celtde-#. (4.4) 
where c and d are complex and |c| > |d|, then 
z;(C;) = (e+y,d)f,4 (d+-y;é) C; (4.5) 


: ; = : i aa ye : 
is the required mapping function, having no singularities in |f;| > 1. By 
using these mapping functions and writing 


OM (z;) = f9(C;), QB’ (z;) = FFE,), (4.6) 

the equations (3.15)-(3.18) may be written 
SF BDAPSMUOASUO+7f40H =O on C,, (4.7) 
ASM) + MSO +ASMO+ nef) = Mit) on C,, (4.8) 
SIO+ASIOHIO+7ALFD = 9 on C, (4.9) 


AS O)+mSFO+ASEO+ue SFO = —M(t) on, (4.10) 
where M(t) = M(t,?): equations (3.10) and (3.11) may be written in a 
similar way. 

The complex potentials OF (z;) which correspond to an isolated normal 
load at an interior point of a plate with fully clamped boundary, may be 
obtained by a method similar to that used by Green and Zerna (3) in the 
corresponding problems in plane strain or generalized plane stress. Thus 
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in case (i), if the force is of magnitude Q and acts at the point z = c, 


and if ¢, = c+y,é = 2,(e;) = (I+y,e; (4.11) 


then 


are) = (1) 
=e \ 


v (*; 07) Haley C4) —2}* log(Z, ~iy)| +Ay{2,(01)—¢,}* log(f,—«,), (4.12) 


Y2 


(24) Ayz,(L1) —€,}* log(f, —é,)+ 
+Y/ 


and Q3(¢,) is obtained by interchanging suffixes 1 and 2 in (4.12). The 
constants H, and H, satisfy (3.20) and the equation 
Q = Sewn _ a9 47_ a7). (4.13) 
3 ln Ye J 
These equations, and the corresponding ones for plates with clamped 
circular and elliptic holes, are derived in (4). 


5. The mixed boundary-value problem for the half-plane 
Consideration is here given to equations (4.7) and (4.8), in which C is 
now the real axis and the functions f9(¢) are holomorphic in S*+ (the upper 
half plane). 
Define two functions ®(Z), ‘Y(¢) sectionally holomorphic throughout the 
whole plane including the point at infinity, by the relations 


SUD)ASFAUE), { in S* 
® => i > (5.1) 
” —AFO)—7eFX0), ¢ in S- 
, vil()+rvefXC), ¢ in S* 
¥(0) = p (5.2 
" | —FE—-fXt), fin S- . 
then OZ) = —¥ (2), (5.3) 
and (4.7) shows that 
+(t)—O-(t) = 0 = ¥+()—¥-(t)_ on G,. (5.4) 


In order to satisfy (4.8), the left-hand side of that equation is rewritten 
using ®(£) and ‘¥({). The expressions for f{({) and f}(¢) are 


(viva) f 20) = ¥(L)—veO(C), ¢ in S*, (5.5) 
(ve—vd S20) = ¥(L)—yi OCC), in S*. (5.6) 

By using these equations and the relations obtained from (5.3), namely, 
oH =—-¥-(), OH)=—V¥+) ond, (5.7) 
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equation (4.8) becomes 


—_ Aa ¥+()—y20+(0) AF) —y, ©*()} + 
Y1—Y2 


l . — . a 
4+ ——— {p, (7, ‘¥-(t)—-®~(t)) —p,(7, ¥-()—®-(t))} = M(t) on Cy. (5.8) 
7 Fe 
Equation (2.3) shows that 


= 4 re. 2a, Vi¥2_ 4s 
Atha = ——, == -—, 
Y%1 Ye a3 Y1¥2 4% 
and these relations combined with the definitions of A; and yp; in (3.8) 
enable (5.8) to be written 


A{@+(t)4+0-(t)} + BY+(t)+0¥-(t) = M(t) on G,, (5.10) 


where the complex constant A and the real constants B and C are given by 


(5.9) 


> 
“I 





(yi; +ye) , as " 
A a, +4, B= a,- —, C = a3y1¥2—4,. (5.11) 
Y1 72 Y1 V2 


Equation (5.7) shows that the complex conjugate of (5.10) is 


A{¥-(t)+'¥+(t)}+ BO-(t)+-Co+(t) = —M(t) on G,. (5.12) 
In order to solve (5.10) and (5.12), numbers v are sought which make 
the equation [(5.10)-+- (5.12)] have the form 
KO+(t)+0-(t) = U(t) onc, (5.13) 
with O(f) = O(f)+ OF (Z). (5.14) 
It is easily seen that vy must satisfy the equation 
BivA C+vA © 
A+vC A-+vB 


which on eliminating v gives two values @,, 6, for 0, namely, the roots of 


6, (5.15) 


A@—(B+C)O+A = 0. (5.16) 
Equations (5.13) and (5.15) show that the corresponding values of v and 
Kk are 
] Ab—C 
a ; == 5.17 
,. -” on — 
and that « satisfies the equation 
(AA — BC)x?+-(B®?+ C2?—24A)x+(AA—BC) = 0. (5.18) 


By using (5.9) and the conditions |y,| < 1, |y,| < 1 it may easily be 
shown that (B+ C)* > 444, which ensures that @,, 6, are different and 
Ky, Kg are real. When the material is orthotropic then not only are the 
roots «,, K, real, but also positive; the proof depends on the conditions 
imposed on the elastic constants by the positive definiteness of the elastic 
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energy function, and for orthotropy these conditions are (see Gladwell (4)) 
a, >a,> 0, a, > 0, a,(a,+a,) > 2a}. (5.19) 
The author believes, but has been unable to prove, that «,, x, are positive 
in the general case. This certainly is true for the numerical values derived 
from the elastic constants given for ethylene diamine tartrate and lithium 
sulphate, which are the only two substances known to the author which 
have the most general type of anisotropy (4): 
Returning to (5.13) it is seen that, in principle, the problem is now 
solved; the two roots 6,, 6, of (5.16) give two functions 0,(2) which satisfy 
Q,() = O(0) +6, ¥(0), (5.20) 
QF (t)—O;() =0 onG, (5.21) 
«;Q7 (t)+0; (t) = U,(t) on C,, (5.22) 
Ut) = (6; M(t)+ M()/(A0,— B), (5.23) 
and these equations may now be solved by known means to give ®({), 
(0), FU(G,), and f2(¢,). In particular, when M(t) is given by (4.10), it may 
be shown that 
Ut) = —(e + ISTO+SEO+O(STO+ ref FO)}- (5.24) 
In order that the stresses should be of an order of infinity of degree less 
than 1 at the end points of C,, it is necessary that ©,(¢) should be bounded 
there. Restricting attention to the case when the roots x,, x, are positive 
and writing X,({) for the fundamental solution of the equation 
x, X#(t)4+X;(t) =0 ond, (5.25) 


which is bounded at the end points of C, and of lowest degree there, the 
general solution of (5.22) which is bounded both at the end points of C, 


and at infinity is ame XA) f U,(t) dt | 
: J XF(HE—-O) 
Using (5.25), and writing [, for the contour Cj +-C,, equation (5.26) may 


be written @40) — X,(2) i,(t) dt 
0 Bile N) J XKOHO 





(5.26) 


» > 
aTTUK 
1 





(5.27) 


6. The mixed boundary condition for the infinite plate with a 
circular or elliptic hole 
The problems for plates occupying these regions may be solved in a 
way exactly similar to that of the previous section, the only difference 
being in the definition of the functions ®(f) and ‘Y(¢). If the circle C in 
the ¢-plane has radius a, then ®(f), ‘Y({) are sectionally holomorphic 
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functions defined throughout the {-plane, including the point at infinity, 
by the relations 


M(Z) — (fC) +F2(0), fin S+ ( C| > a) 

7, Fe0)—7eFU@D), Fin S- (\f| < 

W(t) = (WiS(l)+yef 2(2), Cin S+ | 
\—F9(a2/t)—F9(a2/t), Cin S-J 


The equation (5.3) is replaced by 


(6 
i (6.1) 


a) 


(6.2) 


M(a?/f) = —‘¥(f), (6.3) 


but the remainder of section 5 applies unchanged. 


7. Constant bending moment and shear on part of the real axis 


The problem to be solved now, as an application of the theory given in 
section 5, is that of a plate occupying the upper half-plane, clamped along 
part of the real axis and acted on by constant bending moment and shear 
along the remainder (this is type (b) of the Introduction) which it is 
assumed lies wholly in the finite part of the plane. (Note that if this is 
not assumed then at infinity there are higher infinities in the displacements 
than normal (see section 3).) 

If the bending moment M and shear S are constant then (3.4), (3.5) 
show that K+izF = Mt—HiSt+iOt+e, (7.1) 
and therefore M(t) may be taken to be 

M(t) = 3(Mt—JiSt®)/8h?2. (7.2) 


For simplicity C, is taken to be the single segment |t| < a. If x,, x, are 
positive, then the fundamental solution X,(f) of (5.25) is 


Xj(f) = (C/a—1)*+#94(Z/a+-1)§-*5, a; = (logx;)/27, (7.3) 


where that branch is taken which is holomorphic in the plane cut along 
C,, and at infinity has the form 





X,(f) = C/a+-O(1). (7.4) 
It may easily be shown by contour integration that 
X ;(¢) M(t) dt 
2ni_ J X,()(t—Z) 
I; 
3 


= —jpjal(2ME—iSC*) +X, (C)[iSal—2(Ma+ Sata,)}} (7.5) 
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and from this result 0,(¢) may be obtained by using equations (5.23) and 
(5.27). After simplification it is found that 
3 


(-) { SF Oreo 
(3) = — TeRA(0,—0,) * 


<{2M(1+-0,)[f—aX(0)]+48(1—0,)[22@—aX ()(L+-2iaay)}}, (7.6) 
with a similar equation, with the suffixes | and 2 interchanged, for 0,(Z). 
Since equations (5.5), (5.6), and (5.14) give 


fiGr) = Q4(%), Fale) = Q4(z2) (7.7) 
as linear combinations of ©,(f,), @,(¢,), O2(¢,), and ©,(¢,), the problem is 
now solved. 

Problems concerning isolated normal loads may be solved in a similar 
way; the contour integrals being evaluated by methods similar to those 
used in (1) for isotropic plates: details may be found in (4), where the 
uniqueness of the solutions is also considered. 
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SHRINK-FIT STRESSES BETWEEN TUBES HAVING 
A FINITE INTERVAL OF CONTACT 
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(Dept. of Civil Engineering, University of Bristol) 


[Received 28 November 1957. Revise received 15 May 1958] 


SUMMARY 


The eveluation of shrinkage stresses when the contact interval between two tubes 
is infinite has already been discussed elsewhere (1). Part I of this paper discusses 
the evaluation of shrinkage stresses when the contact interval between two infinite 
tubes is finite. The method considers each tube separately and uses relaxation 
methods to find the radial displacement distribution along its length due to a unit, 
axially-symmetric, radial pressure acting on a small length of each of the two 
surfaces which are in contact. An integral equation is then solved embodying these 
two displacements to find the shrinkage-stresses. Part II discusses the problem 
when both tubes are finite (e.g. collar on a shaft) by solving the elastic equations 
for the stress-functions ¢ and y% (2) with special conditions at the contact interface. 


Part I 


1. We suppose that a tube, shown in section in Fig. 1, occupies the region 
of space defined by 6 <r <c¢ and —w < z < w and is subject to a unit 
axially-symmetric, radial pressure on a small length of its inner surface. 
The radial (uw) and longitudinal (v) displacements must then satisfy 


eu lew u 1—2o0 d*u 1 dv 
(1—o(< — |+ “om - at > == 0 


2 d22 ° 2 Ora 


at ts or 
(1) 
and 09 ae gy) er, oh eens Pa — 


ey 1—20/e*v lév\ 1 du. 1 eu 
Or2 ' r Or] ' 2 0réz' 2rdz 


subject to the stress boundary conditions; these may be written in terms 
of displacements according to 


l+on ou o ou uu. ov 
— 7 = —+ -{|— nee - 
E Or 1—2e\0r 'r § @ ‘ 
; ; (2) 
l+on l/eu aw 
and — 27 om —|-— 4 - 
E 2\ez Or, 


where o is Poisson’s ratio. 

This system of equations and boundary conditions may be solved by 
the relaxation method (2) in which the boundary DC is moved to the 
right as computation proceeds so that residuals upon it are of no conse- 
quence within the accuracy required. The radial displacement distribution 
on r = 6 will be denoted by U. If a second tube, occupying the region of 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 1, 1959) 
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space a <r <b and —w < z < © is subject to a similar unit pressure 
acting on a small length of r = 6 in the negative r-direction, the displace- 
ment distributions may be found in a similar way. The radial displace- 
ment on r = 6 for this inner tube will be denoted by U’. 

Suppose now that the outer radius of the inner tube is not 6 but 6+-3(z), 


— 
y 





—_ 
st-4--b 7] - 


+ 


















































where 8(z) < 6 and is zero unless —l < z < /; this condition defines the 
contact interval. Then, if the outer and inner tubes are assumed to be 
fitted together, the shrinkage-stress p(z) is given by 


I 


8(z) = | [Ue—€)—U"(e—€)]olE) dé, (3) 


“1 

where U and U’ have already been found, and 4(z) is specified. The 
integral in (3) can be replaced by a summation formula and, if the range 
of integration is divided into n equal parts, (n+ 1) simultaneous algebraic 
equations are obtained each involving all the (n+-1) shrinkage-stress 
values. Such a system is easily solved on an electronic computer although, 
if n is less than about ten, hand computation is feasible. 


2. In the numerical example which follows a tube of internal radius 
one unit and external radius two units is shrink-fitted to a solid shaft of 
unit radius and of the same material (o = }), and a contact interval of 
—1<2z<1 with the ‘lack-of-fit’ 5(z) taken as a constant such that 
ES = 10°. With mesh-lengths of 0-1 and 0-2 in the z- and r-directions 
respectively and the unit pressure applied over a length 0-1 the distribu- 
tions of 10°£U and 10°EU’ are given above the line in Figs. 2 and 3, 
respectively. The accuracy of these distributions may be checked in the 
following way. Firstly, if the unit load were to be applied to the whole 
of the inner surface of the tube, instead of a small part of it, the radial 
displacement of any point on this surface would be given by an integra- 
tion of the values given above the line in Fig. 2; the value thereby obtained 
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is 1,898. But, if the unit pressure is applied over the whole surface it can 
easily be shown (by putting V = 0 and wu independent of z in (1)), that, 
for plane strain, the radial displacement (x 10°Z) at any point of the 
inner surface is 1,875. There is, therefore, a small finite-difference error 
in the relaxation solution. A similar treatment of the shaft gives 625 as 
the exact value of the radial displacement (x 10°) on the outer surface 
compared with 646 given by the relaxation solution. 

To find the shrinkage stresses the integral in (3) was replaced by a 
summation based on the trapezoidal rule with a mesh-length of 0-1; the 
eleven simultaneous algebraic equations were then solved for the eleven 
unknown p-values. These are presented above the line in Fig. 4. If the 
contact interval had been infinite and, therefore, U, U’, and p independent 
of z, it is easily shown from (3) that the uniform shrinkage pressure would 
have been 400. To refine the accuracy of the evaluations of shrinkage 
stresses it was first necessary to obtain more detailed radial displacement 
distributions and this was done with mesh-sizes 0-05 and 0-1 in the z- and 
r-directions respectively with the same radial load as before now acting 
over a length 0-05. By the principle of Saint-Venant this redistribution 
of load will not affect the displacement values at a sufficient distance 
away and, in the relaxation computation, the finer mesh was extended as 
far as was necessary. The changed displacement values are given below 


the line in Figs. 2 and 3. In recalculating the shrinkage stresses close to 
the end-point of contact it was assumed, because they were close to the 
value obtained for infinite contact, that the values previously obtained 
on the coarse mesh for —0-5 < z < 0-5 would not be significantly changed. 
The finer mesh values for the shrinkage stress close to the end-point of 
contact are given below the line in Fig. 4. 


3. In the foregoing method the greatest effort is expended in obtaining 
the radial displacement distributions and it can be noted that these dis- 
placements are more easily obtained when the ratio of tube thickness to 
radius is small. For it can be shown (3) that if a thin-walled cylindrical 
tube, of thickness ¢, is subject to an axially-symmetric band of pressure, 
the resulting displacements can be obtained by discussing the bending of 
a beam of unit width, depth ¢ and relevant second moment of area 
t8/R(1—o*) on an elastic foundation whose reaction constant is EHt/R?, 
where # is the radius of the middle surface of the tube. Tables of functions 
are available (3) for use in solution of this second problem so that displace- 
ment distributions are readily obtained. Although this method of obtain- 
ing displacement distributions requires the ratio of tube thickness to 
radius to be small it is interesting to consider its application in violation 
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of this requirement, but for the case of infinite contact between the two 
tubes considered in the example above. Thus the ‘reaction-constants’ for 
outer tube and shaft are {EF and 4£ respectively and, if 5’ and 5” are the 
(uniform) radial displacements of tube and shaft and 5’+-5” = 4, the lack- 
of-fit, then the shrinkage stress p = {#3’ = 45” which gives p = 0-4E6. 
This is the result found previously by a different method and this simple 
way of finding shrinkage stress for infinite contact is valid for tubes of 
any radius. 


Part II 

4. The method described in Part I applies generally only when both 
tubes extend to infinity in the z-direction so that one computation of the 
radial displacement suffices for any position of the load on the contact 
surface. Prohibitive labour is involved in calculating the different dis- 
placement distributions for each position of load application which would 
be necessary if either tube was finite. If the ratio of thickness to radius 
is small such a procedure is feasible on the beam on elastic foundation 
analogy, but the important problem of a collar shrink-fitted to a shaft 
does not come within this category. Here, a method of solution is given 
for this type of problem. 

If an axially symmetric solid of revolution is subject to axially-sym- 
metric loading, the stress-components can be written in terms of two 
stress functions, ¢ and w&, as (2) 


™ l a. ~ 1a 
rr (#, ad — - (4 (l1—o)d), 22: — % 
‘er oy; + r Or 
6 1 le "7 
66 = 2 4 -(f~+(1—o)d) and z= -— 4 
ter r Oz 
The stress functions themselves satisfy the simultaneous equations 
a2 | 4 a2 22 l c a2 a2 
Se hs eB Fs Od (5) 
Gr? or Or * G28 OF rir’ i 2 
. u Ls wT 
Now = €99 = —|00—o(rr+2z)]| 
r E 
where E£ is Young’s modulus, so that, using (4), we have 
l+o 
= (+(1—o)¢). (6) 


Suppose the two tubes (Fig. 5) to be in contact along the whole length of 
the smaller tube (A B) and that the initial lack-of-fit is 5(z). If 27 = 0 on 
AB and s is measured round the boundary we have from (4) that a/és = 0 
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round ABCD and EFGH. Thus y is a constant round these boundaries 
and can conveniently be taken as zero. Hence from (6) 


1—o? 
= (7) 


on the boundaries and, in particular, on the contact surface AB. If u(z) 
and u’(z) are the radial displacements on the contact surface for the outer 


u 
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and inner tubes respectively and ¢ and ¢’ are the corresponding stress 
functions, we have 


u(z)—w'(2) = (2) = =F gg’) (8) 


where R is the rad‘us of the contact surface. 


5. On a square mesh, of size h, the governing equations (5) can be 
expressed approximately as 


$i +$3+(1—$h/ro)be+ (1+ th/ro)p,—4hq = 0 . (9) 
and y+ yYg+(1—gh/ro)ho+ (1+ dh/ro)yhy— Sho + 2bo—$1—$3 = 0 
All external boundaries are stress-free so that (4) and the %-governing 
equation, which is not otherwise satisfied, can be used to eliminate ficti- 
tious stress function values from (9) on these external boundaries (2). 
On the interface, if the first equation of (9) is applied to the outer tube, 
¢, is fictitious; similarly ¢, is fictitious if this equation is applied to the 
inner tube. At the interface 7r = p(z), the unknown shrinkage stress and 
this condition, stated for both tubes in terms of ¢ and ¢, from (4), intro- 
duces in addition a fictitious ys, related to the outer tube, and a fictitious 
ys, related to the inner tube. At the interface, however, it is necessary to 
satisfy the equations (9) for both inner and outer tubes and the equation 
rr = p(z) for each tube. It is therefore possible to eliminate the four 
fictitious stress-function values and the unknown p(z) from these six 
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equations. Using also (8) we obtain 
h h ’ ee 
( — sp) betde)+(I + sa) +h4)— 29+ 


1s ae 2 

4 v5 R "Oe AE, oyh bin tee =0 (10) 
1—o? 2R R 4R? 

for the interface. 


6. For the example which has been considered, the outer tube occupies 
the region of space defined by 1 <r < 2, —1 <z <1 and is shrink- 


370!" 370, 370, 300, 360,68, 360, 375,385, 1, S12, 


¢ 





- 04 > 
Fic. 6. Shrinkage stresses. 


_ 


fitted to a solid shaft of unit radius and bounded by —2 <z < 2. Over 
the contact surface 5, the lack of fit was assumed constant and £6 taken 
to be 1,000; o was taken as 0-25. In a preliminary computation a mesh 
size of 0:2 was used in both r- and z-directions but a finer net of 0-1 in 
both directions was subsequently introduced near the contact interface. 
The stress-function distributions are not given here but the derived 
shrinkage stress values are shown in Fig. 6. It is interesting to note that, 
with the dimensions as above, the shrinkage stress value for infinite con- 
tact given by Southwell (1) is 375 and that, for finite contact, this value 
obtains over the greater part of the contact interval. A more detailed 
solution near the end point of contact could therefore be obtained with- 
out involving the greater part of the contact interval with a resultant 
saving of effort. 
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SUMMARY 
Donnell’s equation for thin-walled circular cylinders is replaced by 
v2 ah. -, Ow = a* 4 
V4(V2+ 1)°w+4K a DY q; 


where w is a non-dimensional form of the radial displacement and gq is the distributed 
radial loading. This equation retains the essential simplicity of the original but, 
unlike Donnell’s equation, the accuracy does not decrease as the wavelength of 
circumferential distortion increases, 

1. Introduction 

ALTHOUGH the precise formulation of the linear theory of elastic shells has 
received considerable attention and is still the subject of controversy, it 
remains apparent that when numerical values are required for a particular 
problem it is often necessary, and indeed sufficient, to resort to relatively 
simple equations in order that the computations should not be too 
laborious. For example, as long ago as 1933, Donnell (1) derived a simple 
set of equations which he used to investigate the stability of thin-walled 
circular cylinders under torsion and obtained an acceptable agreement 
with experiment. 

More recently, these same equations have been employed by Hoff (2, 3) 
and his collaborators for the solution of equilibrium problems of the thin- 
walled circular cylinder, using eigenfunctions which are trigonometric 
either along a generator or around a circumference. However, it is well 
known that as the wavelength of circumferential distortion increases so 
does the error in Donnell’s equations and, with this in mind, Hoff (4) 
later determined the limiting values of the parameters beyond which the 
error becomes intolerable. For this purpose the more accurate, and more 
complex, equations derived by Fliigge (5) were used as the standard for 
comparison. 

The purpose of the present paper is to remove the source of the error 
by improving on Donnell’s approximation and yet retaining the essential 
simplicity. With equilibrium problems in mind, it is then demonstrated 
that the resulting eigenfunctions are quite accurate for all values of the 
parameters. 


(Quart. Journ, Mech. and Applied Math., Vol. XII, Pt. 1, 1959] 
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The application of the new equations to the solution of other types of 
problems is not considered. 


2. Statement of Fliigge’s equations 


Fliigge’s (5) equations for the thin-walled circular cylinder, shown in 
Fig. 1, are used later as the standard for comparison of the approximate 
equations. It is therefore appropriate to begin by listing them. 





Fic. 1 





The equilibrium equations for the element shown in Fig. 2 are, in 
terms of the displacements, 
4 9+ 29s) + 3—vy + —_— ’ 
Viw+ 2w gg+-wt (3 —v)v 29+ U rep—4(1—v) tt 299+ 
3 


2 
4. 12(7) (w—vg—vu,) = af (1) 


U xe +4(1—v)u gg +3(1+v)v—vw + 
1 /h\? 
tg als) {W pr2—4(1 —v)(w .99—U.99)} = 0, (2) 


2 
v gg t+3(l—v)v.+3(1+v)u—we + ae) {4(3—v)w ,.9+ 3(1—v)v.2} = 0, 
(3) 


where the subscripts following a comma indicate differentiation. In these 
equations, A is the thickness and a is the mean radius of the shell, v is 
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Poisson’s ratio, V? is Laplace’s operator 


oe e 
_ 
V  - +35 
and the flexural rigidity is D= an 
sieht ~ 12(1—r)’ 


where £ is Young’s modulus. The distributed radial load is q per unit 
area and the non-dimensional distances and displacements are defined by 
the equations ' , , 

f=—-, U=—, VSe-—-, W=—, (4) 


where x’ is the distance measured along a generator, af that, measured 
around a circumference, and wu’, v’, w’ are the displacements respectively 
along a generator, a circumference, and a radius. 

Fliigge further shows that significant simplifications can be made for 
thin shells where the parameter (h/a)* can be neglected in comparison 
with unity. In fact, the equilibrium equations (1), (2), and (3) can then 
be rewritten in the more convenient form 


V8w+ 4K tw pee t+ 2YW przere + 8W reer t ‘ 
4-2(4—v)w 9909 +-2(2—v)2 2299+ 200 ppo900 + p900 = Ah (5) 


1 [h\? 3(1—v) 1 
Viu = VW 222 — W299 + wala) (1 20000—" zxnee + ( v) ale y Waa 





2 12 
a (6) 
= 249 a0} tae 755) 2eaneat +2 eae) (1) 
In equation (5), the parameter K is defined by 
4K* — 21—w(7). (8) 


To complete this list, the stress resultants in terms of the displacements 
are, according to Fliigge, 





_. 1/(h\? 

N, = alte tey—w)+ ala) xa? (9) 
Ny = j—ilte— wt 75 le y (wpe+)} (10) 
. Eh fa. 1 (h 

Nw = apy \"ot +i5() (wate), (11) 


" Eh |, 1 (h\? 
Me = xen ltt tem ialq) (eae—ualh (12) 
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while the moment resultants are 


D 
M, a te {O27 t+ ¥(W 99+ Ve)+U,}, (13) 
a’ ; j , 
dD. 
M, = ——{wot+w+vw,,}, (14) 
a 
D S 
M4 (l—v){w »+v,}, (15) 
a 
1D 
Mo, = —5 —(1—v)\2wgt+v.,—uUy}, (16) 
2a . , ; 


and the shear stress resultants which accompany these moments are 
D 


Q, , 2 (WO ree TW 700 T Ure ( l —v)U 99+ d(1 t v)v 26} (17) 
as 


Qs 


; {W ggg t+ Wp, + gt (1—v)v,,}. (18) 
3. Donnell’s approximation 
Fliigge’s equations are often avoided because of the laborious calcula- 
tions which are associated with them. On the other hand, the approximate 
equations given by Donnell (1) have met with fairly wide application. 
Donnell’s approximation will not be discussed in detail here. In a sum- 
marized form it is equivalent to the neglect of: 
(i) the terms in (h/a)* appearing in equations (9) to (12) for the stress 
resultants, 
(ii) the terms in wu, v, w and their first derivatives in equations (13) to 
(16) for the moment resultants, 
(iii) the terms in the first and second derivatives of u, v, w in the 
equations (17) and (18) for the shear stress resultants. 
With these modifications, the equations of equilibrium in terms of the 
displacements become 


i as 
V8w+4K4w , pp’ (19) 
Viu = vw pp,—W 299, (20) 
Viv = (2+v)ws29+ © p00- (21) 


4. The improved approximation 

It is well known that as the wavelength of circumferential distortion 
increases so does the error in Donnell’s approximation. However, a close 
examination of Hoff’s results (4) reveals that this error is only significant 
for that particular class of eigenfunctions which becomes increasingly 
independent of position x along a generator. Under such circumstances, 
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it is clearly inadvisable to neglect the terms in v and w and their derivatives 
in the equations (13) to (18) for the moment and shear stress resultants. 

If the above-mentioned terms are retained and the usual substitutions 
made, it is possible to obtain a new set of equilibrium equations in terms 
of the displacements but, unfortunately, they are just as difficult to solve 
as those derived by Fligge. However, in an approximate analysis such 
as this, there is no particular virtue in striving for the elegance of a system 
of equilibrium equations which are necessarily consistent when they are 
expressed in terms of either the displacements or the stress resultants. 
Instead, the aim is for the greatest simplicity at each stage consistent 
with the degree of approximation. With this in mind, the following 
system of equilibrium equations is proposed, 
222 “vA ’ 
Viu = vw p.,—W 299; (23) 
Viv = (2+v)w 9+ W 960> (24) 
the justification of which is given below. 

If the Fliigge equilibrium equation (5) is taken as the standard for 
comparison, it is noticed that the term 


V4(V2+ 1)?9w+ 4K tw 


—2(1 —v)(W prrrer— “270000 — © 2206) (25) 


is neglected in arriving at equation (22). Now, if a further comparison 
is made, and the usual ‘order of magnitude’ arguments are applied, it is 
seen that the constituents of equation (25) are of the same order as those 
neglected in obtaining Donnell’s equation (19). Thus, equation (22) corre- 
sponds at least to the same degree of approximation as Donnell’s equation 
(19). Furthermore, as pointed out at the beginning of this section, the 
error in Donnell’s approximation only grows significantly for the class of 
eigenfunctions which becomes increasingly independent of the position x 
along a generator. Therefore, in arriving at equation (22) account is taken 
of the fact that terms in w», are more important than terms in w,,—hence 
the final neglect of the term (25). 

The remaining equations (23) and (24) are identical with those of 
Donnell and there is no necessity here for their modification. 


5. Eigenfunction solution with trigonometric expressions along a 
generator 
Hoff and his collaborators (2,3) have shown that solutions to many 
equilibrium problems can be conveniently obtained with the aid of eigen- 
functions which are trigonometric either along a generator or around a 
circumference. 
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Following an identical procedure, the solution of the homogeneous 
equation (22) is obtained for the case when the eigenfunctions are trigono- 
metric along a generator, i.e. 

w = eP? cos nz, (26) 
where n is real. When this expression is substituted into equation (22 
one obtains the auxiliary equation 


(p?—n?)?(p?—n?+-1)?+4K4n* = 0 (27) 


which is quartic in p?. From this, the following expressions can be written 
down for the roots p7, p3, p3, and pj. 


Pi, P2, Ps, Pi = n° —44+($+ 21K’)! (28) 


and, since it often occurs that n > 0-17 we have the approximations 
p?, ps, ps, py = n?’—F+ (xn ae +1| Kn — Ss (29) 
ae a 16Kn} 16Kn)| 


which in the limiting case of n 0-17 are as inaccurate as neglecting (h/a)* 
in comparison with unity.t The roots of Donnell’s equation (19) are 

Pip: Pip: Pip Pip = ”?+(Kn+iKn), (30) 
and since Donnell’s equation is rarely valid for n < 0-17 it is seen that 
there is no loss of simplicity here. 

For purposes of comparison, some numerical values of the roots calcu- 
lated from equation (28) are tabulated in Table 1 below along with values 
of the roots obtained from Donnell’s equation and those from the more 
accurate Fliigge equation (taking Poisson’s ratio as vy = 0-3). Following 
Hoff (4), three values are assumed for K, namely 5, 10, and 50. When 
K is smaller than 5, then the assumptions of thin shell theory are invalid 
and (h/a)? cannot be neglected in comparison with unity. When K is 
greater than 50, the shell is so thin walled that it can hardly fulfil struc- 
tural requirements. 


+ The limiting factor is derived from.the consideration that 


(4-4 2iK*n2)t Kni(l- ‘(i ; _ 


, l = l 
= Kn £)(1+ sipqa—~) F taxi! sraxeet~)): 
where the term 1/512K‘n* can be neglected provided that it is sufficiently small in com- 
parison with unity. On substituting from equation (8), we find that 
ee aw 
512K‘*n* 1536(1—v?)n*’ 
but we have agreed that (h/a)* can be neglected in comparison with unity and so 1/512K‘n* 
ean likewise be neglected in comparison with unity provided that 
1536(1—»*)n* > 1, 
or, for vy = 0-3, provided that n > 0-17. 
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TABLE 1 


Roots when w is trigonometric along a generator (w = e”® cos nz) 





K n Pu Pe Pa Ps 

Present 0°0506-+ 0°04941 0°0025-+ 100008 
Fliigge 0°0510+ 0°04901 0°0025+ I-0000t 
Donnell 0°2458+ 0°10177 O1018+ 0°2455% 
Present 0°5316-+ 0° 36961 o1825+ 1°0767% 
Fligge 0°5330+ 0° 36621 o1812+ 1°07831 
Donnell 0°7824+0°31951 O°3241-+ O°7714t 
Present 2°5441 + 0°98011 1'0520+ 2°3704% 
Fligge 2°5469+ 0°9694i 0455+ 2°3749t 
Donnell 2°6278+0°95141 1:0962+ 228071 
Present 12°3921 + 2°01722 7°7414+ 3°22908 
Fligge 12°3952+ 1°9456t 7°7458+ 3°29428 
Donnell 12°4120+ 2°01421 7°7689+ 3°21801 





Present o1012+ 0'09871 o0100+ 1-00021 
Fligge 0°1014+ 009851 o-o100+ 1°00021 
Donnell 0°3476+ 0°14391 0°1440+ 0°34732 
Present 0°9153+0°5133% 0°3619+ 1°29821 
Fliigge O'0157+0°51218 o°3611+ 129871 
Donnell 1°1026+ 0°4535% 0°4567+ 109481 
Present 3°5360+ 1°41318 1°4640+ 3°41311 
Fligge 3°5371 + 1°40971 1°4613+ 3°4146t 
Donnell 3°5963 + 1°39031 1°4923+ 3°3506% 
Present | 14°5372+3°4394? 7°0534+ 7°0887i 
Fligge 14°5388+ 3742091 7°0534+ 710118 
Donnell 14°5535+ 3°4350 7O7II+ 7°O7IIt 


Present 0°5254+0°37418 o'1818+ 108121 
Fligge 0°5254+0°37401 018184 108121 
Donnell 0°7769+0°32171 0°3218+ 077681 
Present 2°3723+ 105121 0-9808+ 2°54251 
Fliigge 2°3723+1°0511t 0°9807+ 2°5426% 
Donnell 2°4585+ 101691 1°0183+ 2°45508 
Present 7°7964+ 3°20651 3°2293+ 7°74148 
Fligge 7°7965 + 3°2063% 3°2290+ 7°7415% 
Donnell 7°8237+3°1954! 3°2409+ 7°71391 
Present 26-2692 + 9°51681 10°9573+ 22°8158% 
Fliigge 26-2695 +9°5158% 10°9566+ 22-8163 
Donnell 26-2776 + 9°51381 10°9616+ 22°8069% 

















It is seen ‘rom Table 1 that for all values of n the numerical results 
obtained from equation (28) are quite close to those from the more 
accurate Fliigge equation (5). The accuracy does deteriorate as K de- 
creases but in connexion with this it should be noted that the value K = 5 
corresponds, for example, to a cylinder which is five feet in diameter with 
a wall thickness of nearly two inches. The error in Donnell’s roots is 
intolerable for the lower values of n. 

The roots of Fliigge’s equation (5) which are given in Table 1 above 
were calculated with the aid of an electronic digital computer but, in the 
absence of such a device they could have been approximated by 


Pr = (l+e)p, (31) 
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where 
(1—v)n?( p*+ p?—n'‘) 
2p*{2(p?—n?+-4)( p?—n?+ 1)( p?—n?)—(1—v)n*(2p*+-1)} 
from Newton’s approximation. In the case when n increases indefinitely 
it is easily checked that the error term ¢ tends to zero. On the other hand, 
when n = 0, both equation (22) and Fliigge’s equation (5) reduce to that 
defining the behaviour of a circular ring. 





6. Eigenfunction solution with trigonometric expressions in the 

circumferential direction 

When the eigenfunctions are trigonometric in the circumferential direc- 
tion, the displacement w is taken to be 

w = e?* cos né, (32) 

where » is an integral number for a closed cylinder. Substitution of this 
into the homogeneous differential equation (25) yields the auxiliary 
-quati ome ‘ , as 
1 eens (p?—n?)?(p?—n?+1)?+4K4p* = 0 (33) 
which is again quartic in p?. The roots p?, p3, p3, pj are found to be 


F _ (2n?—1)\} 





p?, p} = n*—44 iKe| 14+ 


By 
(34) 
Pi, pi= n—4 ine| 1s (1 SCE). 
The roots +pi:p, +Pep, +Psp, +Pap Of Donnell’s equation (19) are 
calculated from the quadratic equations 
Ppt(Uti)Kpp—n? = 0. (35) 


Again, for purposes of comparison, some numerical values of the roots 
calculated from equation (34) are given in Table 2 below along with values 
of the roots obtained from’ Donnell’s equation and those of the more 
accurate Fliigge equation. The same values of K are taken as for Table 1, 
while » assumes the values 0, 1, 2, 3, 4, and 10 appropriate to a closed 
cylinder. 

As before, it is seen from Table 2 that for all values of n the numerical 
results obtained from equation (34) are quite close to those from the more 
accurate Fliigge equation (5). There is a deterioration in accuracy as K 
decreases. The error in Donnell’s roots p, and p, is hardly tolerable for 
n we & s 

Again, the eigenvalues p,,, etc., for the Fliigge equation (5) which are 
given in Table 2 were calculated with the aid of an electronic digital com- 
puter, but they could have been approximated by 


Pr = (1+e)p, (36) 
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TABLE 2 


Roots when w is trigonometric around a circumference (w = e?* cos n6) 





K 


Pi Ps 


Ps Ps 











Present 
Fliigge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fligge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fligge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fligge 
Donnell 
Present 
Fligge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fligge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fliigge 
Donnell 
Present 
Fligge 
Donnell 





4°9502+ 5:05021 
4°9854+ 50154! 
§°0000+ 5:00001 
5°0502+ 495021 
5°0848+ 4°91481 
5019+ 490211 
5°3700+ 4°67641 
5°4026+ 4°63872 
5°4261+ 4763598 
5°9358+ 4°31808 
5°9607+ 4°27801 
5°9931+ 4°2893# 
6°7075+ 3°9845¢ 
6°7354+ 3°93748 
6°7605+ 3°96701 
12°4943+ 3°12491 
12°5167+ 3°03681 
12°5194+ 3°12381 
9°9750+ 10°02501 
9°9926+ 10°00761 
10°0000 + 1000001 
10°0250+ 9°9750% 
10°0425+ 9°9575! 
10°0502+ 9°9503! 
10°1779+ 982817 
10°1950+- 9*81021 
10°2038+ 980421 
10°4413+ 9°5943! 
10°4579+ 9°5761t 
10°4683+ 9g°57182 
10°8243+ 9°2922% 
10°8406+ 9°27321 
10°8522+ 9°27198 
15°2678+ 7°4347% 
15°2808+ 7-40801 
15°2909+ 7°4293t 
49°9950 + 50°0050% 
49°9985 + 50-001 5t 
50-0000 + 5000001 
50°0050+ 49°9950% 
§0°0085 + 49°9915t 
50°0100-+ 49°99001 
§0°0350+ 49°96501 
§0°0385 + 49°961 51 
50°0400 + 49°96001 
50°0851 + 4991518 
50°0886+ 49°91161 
50°0902 + 49°91021 
50°1555+49°8455t 
50°1590+ 49°84201 
50°1605 + 49°8405t 
51°0140+ 49°02551 
51°0174+ 49°02201 
51-0192 + 49°0208% 





° 
01019 +0°09798 
0°3669+0°3195% 
0° 3646+ 0°32208 
0°4261+0°36418 
0°9348 + 0-68001 
09509 + 0°65821 
0°9931 +0°71078 
1°7072+ 1°OI4lt 
I*7OOI + 1°02608 
1°7605 + 1°03308 
7°4947 + 1°8745¢ 
7°4812+ 1°92641 
7°5194+ 187628 


0°0502+0°04971 
o1761+0°17018 
0°1758+ 017041 
0°2038+0°19588 
0°4406 + 0°40491 
0°4430+ 0°40238 
0°4683 + 0°42828 
0°8240+0°70748 
0°8227 +.0°7088 
0°8522+0°72818 
5°2678+ 2°56528 
§°2633+2°5744t 
52909 + 2°5707! 


0°0100-+ 0°01008 
0°0347 + 0°03468 
0°0347 + 0°03468 
0°0400 + 0°04008 
0°0850+0°08471 
0°0850+ 0°08478 
0-0902 + 0-088 
O1554+0°15441 
O1554+15441 
0°1605+0°1595% 
1-0140+ 097441 
1°0139+0°9745t 
10192 +0°9792% 














98 L. 8. D. MORLEY 


where 
= (1—v)(p*—n*+n?) 
ba 2{2( p?—n?+-4)( p?—n?+ 1)(p?—n’*) +4K4p?—(1—v)(3p*—n*+ n?)} 





7. Solution of the inhomogeneous equation 

As pointed out by Hoff (4), when the distributed radial load q is repre- 
sented by a double-trigonometric series then the displacements can be 
expressed in the same form. The procedure is so simple that there is little 
reason to make use of approximate equations in preference to that of 
Fliigge, but Hoff goes on to show that the inhomogeneous Donnell equa- 
tion (19) should not, in any case, be used when the cylinder length is 


TABLE 3 


Double-trigonometric solution of the inhomogeneous equations 











K n s AplA AplAp 
5 I ° - 
5 2 ° 1-0000 1°7778 
5 3 ° 1-0000 1°2656 
5 I o°5 ©°9999 10148 
5 2 os 1°O1S7 14520 
5 3 os 1*0020 1°2557 
5 I 10 0°9999 10042 
5 2 1°0 1°0053 1°0833 
iS 3 1-0 1°0095 I-1904 
5 10 1° I‘OO0I 1°0202 
10 I ° - 
10 2 ° 1*0000 1°7778 
10 3 ° 1*0000 1°2656 
10 I os 10000 1°0010 
10 2 os 1°0020 1°0561 
10 3 o's 10030 1°1835 
10 I 1° 10000 1°0003 
10 2 1° 1°0004 1:0060 
10 3 10 10020 1°O417 
10 10 1-0 I°OOOI 10202 














greater than the diameter. However, this restriction does not apply to 
the equation (22) for its solution is nearly the same as that of the more 
exact Fliigge equation (5). If the radial load q is assumed to be 

q = cos sx cos n0 (37) 


then the displacement w can be expressed as 
w = Acossxcosn6, (38) 


where A, s, and 7 are all constants. Substitution of these two equations 
into the Fliigge equation (5) yields 

Pa: (a*/D)(s?+-n?)? 
F (s2-+-n®)?(s2+-n?— 1)?-+ 4K 4s4+ 2(1 —v)(s®—s2n*+-82n?)’ 





(39) 





AN IMPROVEMENT ON DONNELL’S APPROXIMATION 


and when substituted into Donnell’s equation (19) we find 
(a*/D)(s?+-n?)? 





PD (s?+-n®)t+4K4st’ 


and, finally, when substituted into the equation (22) we obtain 


re (a*/D)(s?+-n*)? (41) 
~ (s?-+-n®)(s?-+ n?@— 1)? + 4K 4st 





For purposes of comparison some numerical values of the ratios A,;/A 
and A,.A,, have been computed and are tabulated in Table 3 above. 

For the particular numerical results given in Table 3 it is seen that the 
greatest error is only 1-6 per cent when the inhomogeneous form of 
equation (22) is used in the place of the more accurate Fliigge equation 
(5). The inhomogeneous Donnell equation should only be used with the 
greatest caution for equilibrium problems. 
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SUMMARY 

A convenient new procedure for direct ‘hyperosculatory’ interpolation for 
f(x) = f(v9+-ph) = f, when values are given for f(2;) = fj, f‘(x;) = fjandf"(x) =f; 
at nm equally spaced points x; = a)+7h, 7% {[4(n—1)] to [$n], is derived from 
Hermite’s (3n—1)th degree osculatory interpolation formula, for n 2(1)7 (i.e. up 
to 20th degree accuracy). Certain fixed auxiliary quantities a;, b; and c¢;, which 
are independent of p, f;, f; and f7, are tabulated exactly. The method is an exten- 
sion of the author’s earlier adaptation of the simple osculatory interpolation formula, 
employing both a decomposition and uniqueness property of Hermite’s general 
formula. The remainder term indicates the vast increase in accuracy and permissible 
size of h, when compared with simple osculatory interpolation formulae which are 
the next most accurate. For inverse interpolation the coefficients of the first ten 
powers of r = (f—f,)/hf;, for n 2(1)5, (i.e. up to 14th degree accuracy in the 
direct function) are given in terms of f;, f; and f7. Hyperosculatory interpolation 
is specially suitable in (1) practical problems in astronautics involving rocket or 
missile flight, where the acceleration data, or fj, are available as well as position and 
velocity data, or f; and f’, and where (2) higher mathematical functions that are 
tabulated with their first derivatives, are solutions of simple second order differen- 
tial equations, so that f/ is readily obtained. 


1. Introduction 

IN a series of previous articles (1)—(6), the writer gave formulae and tables 
for osculatory interpolation (direct and inverse, for real and complex 
arguments) which employs both the first derivative as well as the function 
itself at the interpolation points. The advantages over the ordinary type 
formulae that employ the functional values only are in the much greater 
accuracy when the interval of tabulation is the same, and in that, when no 
increase in accuracy is required, one can interpolate over much larger 
intervals, thus saving a large amount of required tabulation of the func- 
tion. Fuller discussion of these advantages is contained in (1), p. 212. But 
every advantage in osculatory interpolation is enhanced to a still greater 
degree when one can utilize the second derivative as well as the function 
and its first derivative at the fixed points. Such formulae will be designated 
as ‘hyperosculatory’. The present paper will be concerned with giving 
(1) a convenient and concise scheme for direct hyperosculatory inter- 
polation for a real function f(x) that is tabulated at equal intervals, with 
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its first and second derivatives either explicitly tabulated or readily 
available, and (II) formulae for inverse interpolation under the same 
conditions. 


2. Special applications 

At this point it is worth considering some of the reasons for giving 
special attention to formulae that employ second derivatives. One can go 
on giving new formulae employing derivatives of higher order beyond the 
second; but there is a reasonable point at which to stop because in practice 
either higher order derivatives are unavailable or troublesome to calculate, 
or too much time and effort on the preparation of papers on interpolation 
yield diminishing returns in proportion to the need. However, there are 
two very special reasons for developing these hyperosculatory interpola- 
tion formulae, apart from the extreme accuracy over comparatively large 
intervals of interpolation mentioned above:t 

(1) In practical problems of astronautics involving rocket or missile 
flight, where interpolation is desired over the largest possible interval h, 
one often has the acceleration which has been determined by physical 
means (e.g. in inertial guidance), in addition to the velocity and position. 
Thus one may take the fullest advantage of the available information by 
using f, f’ and f” in any interpolation or extrapolation problems on 
trajectories, orbits, or any functions that are tabulated with their first and 
second derivatives. This includes many functions arising in guided missile 
work that are given with their first derivatives, and where the functions 
are solutions of various types of second order differential equations for 
which the second derivative is given fairly explicitly, so that f” is almost 
as accessible as f and f’. 

(2) Many higher mathematical functions of widespread importance 
which are tabulated with their first derivatives satisfy rather simple 
looking second order differential equations. Thus it would be almost 
a pity to employ ordinary or osculatory interpolation, when f” for 
hyperosculatory interpolation is so readily available from the differential 


+ To get a crude idea as to how much better hyperosculatory interpolation is by com- 
parison with osculatory or even ordinary interpolation, consider just the simplest 2-point 
formula using f, f’ and f” at each end-point. Suppose that an ordinary six-point formula 
at intervals of h is required for solving a problem to the desired accuracy. An osculatory 
formula using f and f’ at three points would be a confluent case of the ordinary formula 
where the six entries over a total range of 5h could be regarded as coalescing in pairs, into 
three entries of f and f’ at intervals of 24h, so that the interval of interpolation is now 2} 
times as great. But suppose in that same six-point formula there is coalescence by threes 
into two entries of f, f’ and f” at the two end-points, 54 apart. Thus we now have an 
interval of interpolation that is five times as great as in ordinary interpolation and twice 
as great as in osculatory interpolation. 
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equation f"+f'/x+f = 0 in a table of J,(x) with J,(x) = —J%(x) (7), or 
of Y,(x) with Y,(z) = —Y (x) (8), or of any linear combination of those 
two. Also any type of Hankel function that is tabulated together with 
its derivative (9), is well suited to this hyperosculatory interpolation since 
it satisfies a simple Bessel-type second order differential equation. Like- 
wise for interpolation in any kind of Hermite, Legendre, Jacobi, Mathieu, 
or hypergeometric function, where the first derivative is given with the 
function, one should use f” too, as it is so easily obtained from the 
differential equation. 


3. Direct interpolation 

In the case of direct interpolation we wish to find f(x), given f(x;) = f;, 
f'(x,) =fj and f"(x;) = fj at equally-spaced points 2; = 2,+1th, where 
h is the tabular interval, x = x)+>ph, f(x) = f(z»+-ph) =f, =f. For the 
n points x,, i runs from —[4(n—1)] to [4n], the symbol |[m] denoting the 
largest integer not exceeding m. Then the hyperosculatory confluent form 
of Lagrange’s interpolation formula, or in other words, Hermite’s formula 
as far as the second derivative, which gives the polynomial of degree 3n— 1 
having prescribed f,, f; and f/ at « = 2,, may be shown to be given by 
the following formula in terms of the variable p: 


f(xo+ph) > {LE (P)PP)L:AAG (i) fi + VHD) f+ Ban(P), 


i {4(n—1)] 
(1) 
where 


Fp) = 1-31" (i)(p—i) +(— BL") + O{L""()}*)(p—i), (1) 


Gp) = (p—i{l—3Ly"(i)(p—a}, (1”") 
Hp) 1(p—t)?, (1”) 
Ua) ia 
Li""(p) = i ((p—Jj)/(i—J)}, (2) 
J=—lLi(n )] 
and 
sngiany ey) EP 13 Janis 
Rs,,(p) h ei (p—J); / (3m); (3) 
j {4(n—1)] 


where 2_((,-») S € < Typ) in the case of interpolation only, i.e. for 
Lin SX < yp), and € is in the interval (x, x_:(,_1),---, %n)) for any 
location of x, thus covering extrapolation also. Whenever f®”"(€) is difficult 
to estimate, we may approximate A®"f®@"(£) in (3) by A%”, or the (3n)th 
difference of the entries f;. 

The derivation of (1) is straightforward from the general existence and 
uniqueness theorem on Lagrangian interpolation, including all confluent 
forms (10), and is suggested by the appearance of the explicit formula for 
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ordinary osculatory interpolation (10), (11). Thus in any (3n—1)th degree 
polynomial approximation to f(z,+ph) of the form 
[én] ° . 
PP) = > Mart bdp—d)+edp—O fit 
+h(aj,+-b;(p—i)+¢(p—i)?) fj +h*(aj+-b}(p—t) +¢(p—i)?*) fi}, 
the factor of {L\"(p)}* makes the coefficient of f,, fj or fj vanish in P(j), 
P’(j) or P"(j), for 7 Ai. The condition P(t) = f;, leads to a; = 1, 
a, = a" = 0; the condition 


t 











d 3s. oe 
dc’ '?) Feo | i Pers 
leads to 6, = —3L% (i), b; = 1, bf = 0; and the condition 
d? P"(i) ” 
dx? P(p) —_ h2 ~ fi ’ 





leads to a 
c= EL" OLE", = —3L"),_ of = 4. 
This completes the derivation of (1), except for (3), the expression for the 
remainder term R;,,(p), the derivation of which may be found in (10). 
Some idea as to the extreme accuracy of (1) may be had by a glance 
at the value of R,,(p) in say the 3-, 5- and 7-point cases for p = }. 
Then R,,,(4) would be 


—(—1))(4—0)(4—1)}8 
~ (== YF op, 








{(4—(—2))8—(—) 4-94 —D4—2)? ars 
15! 
and (4-(—3))(4—(— 2) -(— 4 —DG—24—3)F nn 
21! : 
or ~ —1-45.10-7A%, 2-13.10-"2A™ and —3-65.10-'"A*! respectively. 
The corresponding remainder terms R,,(4) in ordinary osculatory inter- 
polation turn out to be 


—(—1))(}—0)(4—1)}? 
_ fiend Me d=) ae 








{(4—(—2))(4—(— 1))(4—0)(4—1)(4—2)}? a 10 
10! 


(3—(—3))(§—(— 2))(—(— OB — DE — 204 — 3)? ae 





and 





14! 
or ~ 1-95.10-4A®, 5-45.10-7A™ and 1-74.10-°A™ respectively, which, 
although still impressively small, are not anywhere near the smallness of 
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the R;,,(4). For final comparison, we note the corresponding remainder 
terms R,(}) in plain ordinary Lagrangian interpolation, namely, 
~ GH(-)G-94=)) 43 
3! 


(4—(—2))\(4-(— 4-14 —DG—2) 4 








(—(—3))(4—(—2))(4—(— 1) 4—9)8— DG—2)(8—3) 





ir | 


‘. 

or ~ —},A%, 3,A° and —;%,A’, which are considerably large by com- 
parison. Another idea as to the power of the hyperosculatory interpola- 
tion formula (1) may be furnished from the value of the remainder in 
extrapolation over one full interval h, which is obtained by setting 
p =[4n]+1 in (3). Whereas ordinary 3-, 5- and 7-point Lagrangian 
extrapolation gives remainders ~ A*, A® and A’ respectively, and the 
corresponding osculatory extrapolation gives the considerably smaller 
approximations of jA®, ,,A! and 2-91.10-4A' respectively, here the 
R;,({4n]+1) for hyperosculatory extrapolation gives the exceedingly 
small remainders of approximately 5-95.10-4A®, 1-32.10-®A and 
2-51.10-*A*! respectively.t 

Formula (1) as it stands is very cumbersome, involving large amounts 
of computation for the coefficients of f,;, f; and f;. In this present treat- 
ment we give a more concise method of using (1) for any p, employing 
a very small schedule of precomputed constants, a;, 6; and c; which are 
tabulated below. The use of this new method saves a great amount of 
labour which would be required in using (1) directly. The principle is 
almost an extension of that expounded in (1), pp. 212-14, the only 
difference being in the extra term involving f/, and it will be recapitulated 


+ This present paper is concerned exclusively with this hyperosculatory or doubly- 
confluent Lagrangian interpolation and the advantages are shown in comparison with just 
ordinary and singly-confluent Lagrangian formulae. Thus the table-maker may still be 
concerned with other aids that are either remotely or not at all related to the ideas here. 
For example, the use of ‘reduced derivatives’ which is synonymous with the Taylor series 
represents a completely confluent Lagrangian formula. There, instead of the function at 
n distinct points, we employ the function and its first (n—1) derivatives at a single point, 
which makes it even more accurate than these present formulae. But in practice those 
higher order derivatives are often either unknown or too difficult to calculate. Some other 
types of interpolation aids like modified differences and coefficients in economized poly- 
nomials (the most popular arising from expansions in terms of Chebyshev polynomials) are 
based upon the entirely different idea of approximating under suitable conditions a poly- 
nomial of higher degree by one of considerably lower degree. Thus, if desired, one may 
start with the explicit formulae for (1) obtainable from s, t, u, v,... in the section on Inverse 
Interpolation and, with a considerable amount of manipulation, approximate them by lower 
degree economized polynomials. 
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only in outline here, for the reader’s convenience. Writing (1) in the form 
Lin] ° 
A? 3A3 L{”"(i) 
i h a Lo 3 .. as + Pa. 
Seat eh) = P)} pogo (es (p—t)? 
3/__ (n)" (4 (n)’(4)\2 
Aj( $l; (i) + 6{L (t)} id 
(p—t) 
A3 3A3L (i) 
h S. 
” oa (p—t) 





oo 








iit gets + Ron) (4) 


(in) 
where Lp) = {i (p—)) (5) 
j=—li(n—])) 
nj od 
and A, = | tr (i—)| . (6) 
j--Gin-1)) 


just as before, by setting f(z,+ph) = 1 and using the uniqueness of re- 
presentation of any arbitrary (3n—1)th degree polynomial, by (4), we 
obtain 
Lin) ons 
A3 3A3L\"(i) 
L™(yy)-3 = i SAL" (+) 
(U\P)} 2 bs “ (G7 (p—1)* ‘a 
AK MPO OLE) oy 
(p—t) 
Thus from (4) and (7), it is apparent that (1) may be expressed concisely 
(omitting remainder R,,,(p), and limits of —[}(n—1)] to [4n] in > under- 
stood) in the form 





— 














L~X (4S APB LAM YA & (8) 
where 
ay = + ta Le, (9) 
e =, om (p—i) 

get a; b; 10 
Pe apt Gay’ sit 
t,o rae (11) 

and the fixed quantities a,, b; and c,; may be defined by 
‘ a; = k(n)A}, (12) 
b, = —3k(n)A? L(”"(i), (13) 
and c, = k(n) A3(— 3 L(t) + Of LI""(a)}?). (14) 


In (12)-(14), k(n) may be any suitably chosen constant of proportionality 
which depends only upon n, and whose function is merely to express a4, 
b, and c, as exact integers instead of exact rational fractions; this simplifies 
the computation in (8)-(11). 





106 H. E. SALZER 


The table below gives the exact values of a;, 6; and c; for n = 2(1)7, 
with corresponding accuracy of degrees 5(3)20. 


Table of a;, b; and c; 








a; b; a) 
2-point a : bo 3 “0 6 
a, I b, —3 4 6 
ay 2 by 9 Cy 24 
3}-point lo 16 0% ° Co 46 
1 2 b, —9 Cy 24 
ay 6 by -33 C4 103 
4-point 49 — bo 243 “o 129 
ay 162 d, 243 a —729 
ly 6 b 3 fe 103 
as 12 bs 75 Cs 260 
Bs — 768 by — 1920 C4 5120 
5-point ag 2592 by ° Ce 9720 
a, — 768 b, 1920 a — 5120 
: 12 by 75 ty 260 
as 300 b_s — 2055 C_s — 7697 
ay 37500 b_, 1 21875 C4 3 34375 
a a 3 00000 by — 3 00000 &o —13 25000 
I a 3, 00000 b, -3 00000 ty 13 25000 
ay —37500 | by 1 21875 | 3 34375 
ly 300 bs 2055 cs 7697 
a_s 600 db, 4410 C3 17549 
as --I 29600 bs —4 98Q60 Cs —14 39424 
a_, 20 25000 by 35 43750 ty 112 21875 
7-point ag — 48 00000 by ° fo —196 00000 
a, 20 25000 b, -35 43750 rm 112 21875 
ly -I 29600 b, 4 98960 Ce —14 39424 
a3 600 bs —4410 fs 17549 























The computation of «,, 8; and y, in (9)-(11) is facilitated by the recur- 
rence scheme of first finding y,; by (11), and then 


B, 27 ty (10’) 
(p—?) 

followed by oy Bites (9’) 
(p—t) 


4. Inverse interpolation 

In n-point inverse hyperosculatory interpolation the problem is to find 
x = %+ph from given values of f(x,+ph), f;, fj and f7, i = —[4(n—1)] 
to [4n]. Here the resulting formulae for p are given by the inversion of 
(1) or the direct (truncated) power series for f in terms of p. Since (1) is 
much more accurate than either the ordinary Lagrange formula or the 
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simple osculatory formula, the resulting inverse hyperosculatory inter- 
polation formulae for p would be expected to be much more efficient and 
enable the user to employ very much larger intervals A for sufficient 
accuracy, surpassing even the permissible size of A in inverse simple 
osculatory interpolation (3). The formulae below apply to all the types 
of problems cited under Special Applications in connexion with direct 
interpolation. All formulae express p as a rapidly convergent series in 
powers of a variable r (a suitable linear function of f), the coefficients of 
r™ being very simple functions (and the same for every n) of fixed quanti- 
ties s, t, u. v,..., Which are given below in terms of f,, f; and f; for n = 2, 
3, 4. and 5. Although the direct hyperosculatory interpolation formulae 
forn — 4and n = 5areof the 11th and 14th degree accuracy, respectively, 
the inversion formula for p that is given below does not go beyond the 
10th degree terms, because in numerical work one rarely goes even that 
far, and in most practical problems the first few terms will be sufficient. 
In other words we simply do not employ the coefficient of p™ in (1) for 
n = 4 nor the coefficients of p™, p'*, p'® and p"™ in (1) for n = 5. For 
every n, we define r = (f—f,)/hf, and s = hf,/2f,, and corresponding to 
n = 2,3, 4 and 5, quantities ¢, u, v, w, x, y, z and Z are defined as follows:t 


n=2 
t = [—10( fo—f,) —h(6fo +41) —$h*(3fo —S7)] /Afo, 
u = [15(fo—f) +A(8fo+ Tf) + $4°(3fo—2f7)] /Afo» 
v = [—6(fo—fi)—3A( fot fi)—F( fof) Aho, 


w= ft=y=o2z=—72= 0, 


= 3 
t= [—35(f1—fi)— ACL + 48fo+ Lf) —A(f" 1 —f7)] /L6hfo, 
= [48( f_1—2fo+i) + 13A(f~1— fi) +4 f — 24 fo +7))/ 16hfo, 
[21( ffi) +h(Of 1+ 24fot Of) +0 fF .—f7))/8hfo, 
[—32(f_.—2fotfr)— Aff) hf — 12fo +f) /8hfo, 
= [—15( ffi) A(T +1 6fo+ Th) hf" —fi))/ M6Af, 
y = [24(f 4 —2fotSd+9h( ff + A(F2 — So +/2)]/ L6hfe, 


= i, 


~ 
= 


~ 
~ 


& 


~ 
~ 


nN 


+ The inverse interpolation formula for p in powers of r works even for n = 1, i.e. given 
only f,, f; and ff besides f, if wesett¢=u>-v=weae2=yuz=2= 0. 
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n=4 
t = [—1136f_,—7452f,+9072f,—484f,— 
—h(312f"_,+5832f,+3240f; — 144f,)— 


—h?(24f” ,+972f7—648f7 + 12f3)]/1296hf,, 
u = [3200f_,—3645f,+-445f,+h(828f" , + 1296f,+324/;— 138 f,)+ 
+h?(60f" ,—1458fp—324f7 + 12/2) ]/1296hf,, 
v = [—1592f_,+21627f/,—21384f, + 1349f,— 
h(246f’_,—9477f,—8262f; + 399f,)— 
—h?(6f"_, —2835f5 + 1782f{—33f3)]/1296hf,, 
w = [—3120f_,—2025f,+ 6480f, — 1335f,— 
h(1107f"_,+3888/5+ 2349f; —414f,)— 
-h2(99f" , —486fo—1053/{+-36f;)]/1296hf,, 
- [3573f_, —20169f,+ 17739f, — 1143f,+ 
LA(1O71f" , —6399f,— 6885; + 333f,)+ 
+h*(81f" ,—2673f5+ 1539f; —27f2)]/1296h/,, 
y = [—120f_,+8505f,—9720f, + 1335f,4 
h(54f’ ,+3888/5+ 3402/,—414f,)-+4 
+-h®(18f" ,+-810fp—1134f {+ 36f2)]/1296Af,, 
z = [—1390f_,+5265f,—4050f, + 175f,— 
h(480f" ,—1215f,—1620f;+45f,)— 
—h?(48f" ,—729fo+ 324f{—3f2)]/1296hf,, 
Z = [688f_,—4131f/,+ 3888/,—445f,+ 
+-h(225f" ,—1296f,—1377f; + 138f;) + 
+h*(21f" ,—486f5 + 405/7— 12f3)]/1296hf,. 


~ 


n=5 
t = [9616( f_,—f,)—499712(f_,—f,)+ 
+-h{2592( f’ »+f,)—172032(f" ,+f;)—622080f¢}4 
+h*{192( f” -—f)—24576( f” , —f7)}]/165888hfe, 
u = [—5504(f_,+f,)+720896( f_,+/,)—1430784/,-+ 


+h{—1392(f* .—f,)+ 221184(f"_,—f;)}+ 
+h — 96(f7 o+f2)+24576(f" ,+-f7)—311040f5}]/165888hf,, 
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v = | —32936(f_.—f.)+882688(f_,—f,)+ 
+h{—8976( "9 +fg)+374784(f* ,+f1)+-902016fo} + 
+-h*{—672(f" »—f7) +67584(f* ,—f)}]/165888hf,, 

@ == [18880( f_.+-f,) — 1392640( f_,+-,)+-2747520f,+ 

+h{4824( f” »—fz) —509952( f"_ ,—f;)}+ 
+h*{336(f" ¢+f2)—67584( f” +f) +451008f5}]/165888hf,, 
aw = [41557(f_,.—f,)—679424(f_,—fi)+ 
+h{11514( f"9+f)—290304( f* , +f;) —635040f9}+ 
+h*{876(f” .—f2) —66048( f” , —f7)}]/165888hfo, 
= [—23876(f_.+/2)+1101824( f_,+f,)—2155896f/,+ 
+h{—6195( f*»—f2)+422400( f" ,—fp)}+ 
+h*{—438(f" .+f2)+66048( f” , +7) —317520f5}]/165888hf,, 
| —23215( f_.—f,) + 256640( f_,—f,)+ 
+h{ —6606( f” +-f3)+ 104064( f”_, +f;) + 225504f5}+- 
+h*{—516(f” ,—fg)+28032(f” ,—f7)}]/165888hf,, 
Z = [13388(f_.+/f.)—416768( f_,+-f,)+ 806760/,+- 
+h{3561(f! »—f,)—160128(f_,—f,)}+ 
+h*{258(f" .+-f¢)—28032(f",+f7)+1 12752f5} ]/165888hf,. 

For every n, the formula for p is the following: 

p = r—rs+r3(2s2—t) + r4(— 583+ 5st— u)+-15(1484*— 218%-+ 32+ 6su—v)+ 
+-r®(— 4285+ 8489t— 28st?— 28s*u + Ttu+ 7sv—w)+ 
+-r7(132s*—330s*t+ 1808%?+ 12089u — 12#8 — 72stu— 36s*v+- 

+4u*+ 8tv-+ 8sw—zx)+ 
-+-r8(— 42987 + 1287s5t—990s9t? —495s4u + 495s"tu + 165st3 + 1658°» — 
— 45t?u—45su?— 90stv — 45s*w + 9uv + 9tw+ 9sx—y) + 
+ r9(143088 —5005s*t + 5005s4t? + 2002s5u — 1430878 — 2860s9tu — 
—715s8'v+ 55t4+ 660st2u + 330s8*u? + 6608*tv + 22083 w — 
—55tu2—55t2v—110suv—110stw—55s%x + 52+ 
+ 10uw+ 10tx+ 10sy—z)+ 
+ 0(— 4862s9+ 1944887t— 24024852 8008s%u + 100108%3+- 
+ 15015s*tu + 3003s5v— 1001 st*— 6006s*t2u — 2002s°u? — 
—4004s%tv — 1001 s4w + 2866 u + 858stu? + 858s%tw+- 
+ 858st? + 858s2uv + 2868*2 — 22u3 — 132twv — 6622» — 
— 66sv?— 132suw— 132str — 6682y + 1 low+ lluxe+ 
; llty+1lsz—z)+.... 
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SUMMARY 


A known nine-point difference analogue of Laplace’s equation is rederived for 
rectangular grids in such a fashion that it is readily established as a ‘ best’ difference 
analogue by elementary methods. 


1. Introduction 


A POPULAR technique for the numerical solution of a Dirichlet problem 
involves the replacement of the Laplace equationt 


Ur t+ Uyy, = 0, (1) 
by the nine-point difference analogue 
— 20g + 4(u, + Ug+ Ug+ Uy) +U,+Ugt+u,+u, = 0, 
where 
Ug = U(X, Yo), Uy = W(Xyot+h, Yo), Uz = U(X, YotA), 
Uz = U(Xyp—h, Yo), Ug = U(X, Yo—h), Us = U(ry+h, yo+h), 
Ug = U(%yo—h, Yo th), Uz = UW(xy—h, yyo—h), Ug = ulrot+h, yo—A), 


and fA is the positive grid constant. (See, for example, reference (1), 
chapters 8 and 10.) Also, it is known that equation (2) is a ‘best’ nine- 
point analogue of equation (1), subject to certain assumptions and defini- 
tions (see reference (2)). 

The primary object here is to determine a ‘best’ nine-point analogue 
of Laplace’s equation for which the grid size h, in the z-direction, need 
not necessarily be the same as the grid size d, in the y-direction. 

The concept of ‘best’, as used here, will be completely defined and 
studied in section 3. Intuitively, a difference approximation is ‘best’ if 
no higher order difference approximation of the same type exists. 


2. Derivation of the difference analogue 

Let the grid sizes in the x and y directions be h and d, respectively. 
The following notation will be used: 

Uy = U(X, Yo), U, = U(X +h, yo), Us = U(X, ¥y+d), 


+ Suffixes denote partial derivatives. 
[Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 1, 1959] 
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Ug = UTo—h, Yo), Ug = U(X, Yo—A), u, = u(ro+h, y,+d), 
Ug = U(Xy—h, Yyt+d), Uz = Ulry—h,yy—d), Uy = ulapth,yo—d) (3) 
She dn inns eee (4) 
P m! nt ox™oy™ 


Since h and d are positive constants, let 


h=pd (h>0,d>0,p> 0). (5) 

Since u is harmonic, it follows from (2), p. 426, that 

Ago = —Ag, As» = —2A,, 

As9 = —4Aj2 A,, = —10A,, 

As, uF 3Ag.s Ago an —Aog 

Agy = Agu A;, = Aj, (6) 

As, Aj; Age 15Ag., 

Ass 6A, As vA15 

As» = $A,, Ao, sie —15Ag.¢ 

Ag, 5Ags | 





For future reference, note that by use of (4), (5), and (6), the Taylor 


series about (x, %) may be written 


Up = 


Uy 


Us ‘ 


tls - 


U- = 


o 


Ago 
Ago +Aj9pd—Ao. p*d?—}4A,, p*d®+ Ag, p*d*+ 3A, ,p*d>— | 
— Ay, p*d*+ O(d") 
Ao +A, 4+ Ao2d?+ Ag3d°+Agyd*+Ay,d>+Ay,d*+ O(a") 
Ayo—Ai ea Ao. p*d?+ 3A, .p*d*+ Ags ptd*—tA, ,p'd>— 
—Ao, p*d® + O(d") 
a= — Ao,d*—Ag,d°+Ay,d*—Ay,d°+ Ay,d°+ O(d") 
Ago +d[pAjo+Aoi]+4*[pA,,+Aoe(1—p*)]+ 
+d®[A, o(p—4}p*)+ A, ,(1—3p?)]+ 
+d*[Ao4(p*—6p?+ 1) + A,3(p—p*)]+ 
+d*[Ao5(5p*— 10p?+ 1)+ A, 4(kp>—2p*+ p)]+ 
+d*[Ao(1—15p?+ 1ip*—p*)+ A, ;(p°—2p*+ p)]+ O(a’) 
= Agg+d[—pAjo+ Ao,)+@*[— pA, +Ao(1—p*)]+ 
+d*[A, o(hp°—p)+Ao3(1—3p*)]+ (7) 
+44 Ao 4(p*—6p?+1)+A,3(p?§—p)]+ 
+d*[Ao,5(5p*—10p*+1)+A,4(2p>—p—}fp')] 
+d%[Ao¢(1—15p?+ 15p*—p*)+ A, 5(ep3§—p— > ]+0@) 
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Uz = Ago +d[—pAyo—Agi]+4*[pAy 1 + Aoo(1—p*)]+ 
+d*[A, o(4p®—p)+ Ao3(3p?—1)]+ 
+d*{ Ao4(p*—6p?+ 1)+ A, 3(p—p*)]+ 
+-d*[Ao,(10p?—5p*—1)+ A, 4(2p*—hp>—p)]+ 
+-d*[Aoe(1—15p?+ lip*—p*)+ A, .(p>—Pp* + p)]+ O(@’) 

Ao +4|pAj9—Ag,]+4*{—pA,1+Aoe(1—p*)]+ 

+d[ A, o(p— }p*) + Aos(3p?— 1)]+ 
+d*{ Aq 4(p*—6p?+ 1)+ A, 3(p?—p)]+ 
+4*[Ao,(10p?—5p*—1)+A,4(p°—2p*+p)]+ 





+-d®[Ao9(1—15p?+ 15p*—p*)+ A, ,(2p?—p>—p)]+ O(a’). | 
Now let L(u) == ¥ amy 


Substitution of (7) into (8), and rearrangement yield 
E(u) = Ag o(a%q+ 04+ ogg + 1g + 45 + 01g + 17 + 04g) 
+-d{A j 9 p(ay —a1g +015 — 1g — 47 + 0%) + 
+ Ag s(og— ag + 45+ 0% — 4p — Og) } + 
+d*{Ao[— a, p? + a— ag p?—a,+(1 —p*)(a5+ a+ 47+ 0%) ]+ 
+Aj1 P(as—a%g+0,—a%)} + 
+d*{ 4 pA, o[ —ay p*+ ag p®+(3—-p*)(a1s—ag— a7 +) ]+ 
+Ags[ag—ag+(1—3p*)(as +14 —a7— ag) ]} + 
+d4{ A, 3(p—p*)(ag—ag+a7—a) + 
+ Ag «lay p*+ ag +05 p* +014 + (p*— 6p? + 1)(a5 +g + 27+) ]} + 
+d*{kpA, gla, p*—ay p*+ (p*— 10p?+-5)(a5—ag— a+) ]+ 
+ Ag slag—a4+ (5p*— 10p? + 1)(a5+ %g—a,— ag) ]} + 
+d{A, 5(p>—Pp*+ p)(a5— a+ a —%) + Agel — a p® +0— 





— a p®+a4+(1—15p?+ 15p*— p*)(a,+ag+07+ a)]}+ O(d")) 
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In order to make L(u), in (9), of as high an order as possible in d, con- 


sider the following system of linear equations 

Agta + Ag+ Agtaytastagta,ptag =O | 

Oy — Ag+ Ag —Ag—Az+ ag = O 

Og— gt a5 -+ cig—A7— Og = O 

— ay p? + Og—as p*+ 04+ (1—p*)(ag+ag+a7+04) = 0 

Ag—Ag+az—ay = O 

—ay p? +a p?+(3—p*)(as—ag—a7+a%) = 0 

Og — ag-+ (1 —3p*)(a5+-1g—a7—a%1g) = O 


at P*+ ag+ ag p*+- a4 + (p*—6p*+ 1)(ag+ag+a7+04) = 0 
I 
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The solution to this system, in terms of a, is 
p?—5\, 1—5p*\_ 
Oy = Ag — ie 4) Qe == Ms = a p+ ’ 


‘lite —_ _ = 1 
er QgerGg=a= “re 





Hence, from (8) 

2. (pe—s 
LU) = % Uq+ io%o| SF J (4a +s) + 
p*+1 


1— 5p? ' 
7 (SO) (Ug+ U4) —gyXo(Us+Ug+Uz+Ug) (11) 
while from (9): L(u) = O(d®). (12) 


Hence 


9 \ 
— 


11. (p*-—5 ,. (1—5p? 
woe tl aT] t+ wa) Home sa (Ug+U4)— 


— yo (astag+a+a5) = O(d). (13) 
Finally, letting a, = —20 and discarding the terms O(d*), (13) yields 
the following difference analogue: 


5—p? 5p?—1 
— 20u,+2 es (u,+Us)+2 ~ — }(Upg+ Uy) +(U,+Ug+U,+Ug) = 0. 
1+p? 1+ p- 
(14) 


One may note immediately that for h = d, ie. for p=1, equation 
(14) reduces to equation (2). 


3. Establishment of equation (14) as a ‘best’ difference analogue 

for Laplace’s equation 

It is now assumed that in equation (5), the only predetermined constant 
is the mesh ratio p, that is, p is a fixed positive constant independent of 
h and d. Hence, by (5), / is a function only of d and h>0 as d—>0. 
This assumption also allows discussions which involve h, or A and d, to 
be replaced by discussions involving d alone. This is readily utilized in 
what follows. 

An allowable difference equation analogue, for present purposes, is 
defined to be an equation of the form 

8 


p32 a; u; = 0, 
where the u; are given by (3) and lima; 4 0 for at least one value of 
a0 


i = 0,1, 2,...,8. (For heuristic considerations concerning this definition, 
see reference (2), p. 428.) If, in the preceding discussion, one constructs 
in place of (13), say (13’) . 


> a; u; = O(d") (13’) 
0 
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and, correspondingly, in place of (14), the analogue 

> a; u; = 0 (14’) 
where (14’) is an allowable analogue, then (14’) is called an approximation 
of order (n—1). 

It will now be shown that (14) is a ‘best’ difference analogue in that 
both the following theorems are valid. 

TuHeoreM |. Jf p = 1, then equation (14) is an allowable 7th order ana- 
logue and there does not exist an allowable 8th order analogue. (For proof, 
see (2), pp. 428-30.) 

THEeoreM 2. If p ~ 1, then equation (14) is an allowable 5th order ana- 
logue and there does not exist an allowable 6th order analogue. 

Proof. The fact that (14) is an allowable 5th order analogue follows 
from equation (13). We proceed, now, using proof by contradiction. 

Suppose that there exists an allowable 6th order analogue. Then let us 
consider in detail equation (9). 

From equations (6), it follows that Ag», Ay», Ags, Ais; Aoa Ara Aes: 
A, ;, Ao, may be considered as independent parameters. Moreover, since 
Ao», Ayo, Ao, are independent of d, it follows that one may consider 
Ao Aor» Aro» Ao,2» A1,2» Ao,» Ai,» Aa» 41, 40,5» 41,5, oe 48 independent 
parameters. Now, by assumption, there exists an allowable 6th order 


analogue, which implies that the right-hand side of equation (9) is of 
order at least 7 ind. Since the A, ;’s which occur in equation (9) are those 


which we are considering to be independent parameters, it must follow 
then that 


tg ++ m+ Ag+ ag +a5+ag+a,+a = €, = Old"), (15.1) 
ty — Ag+ 1g —Ag—Az+ag = €, = Ofd*), (15.2) 
Og — Ag+ 15+ Ag—a1p— Og = €, = Old®), (15.3) 
x, P? + aig— 1g p?+- 1g + (1 —p®)(arg + rg argtarg) = €, = Old5), (15.4) 
Xg— g++ A_—ag = €, = Ofd*), (15.5) 
— ay p+ 0 p"-+(3—p"\(0g—ag—0y +04) = 4 = Of), (16.6) 
Og — ag-+ (1 —3p*)(a5+-1g—07—a¥g) = €, = Old*), (15.7) 
¥y P* + og + ag p* + ay-+(p4— 6p? + 1)(as+ag+ 07+ 04) = €g = Od*), (15.8) 
x4 P*—ag p+ (p*—10p?-+-5)(ag—ag—ay- +4) = & = Of), (15.9) 
ay — ag+ (5p*— 10p?+-1)(a5,+-ag—az—ag) = €y9 = Ofd?), 


(15.10) 
— a, p®+ ay— ag p®+a4+(1—15p?+ 15p*— p®)(a5+ a+ 07+) 


= €, = Old). (15.11) 
Recall, of course, that O(d") represents a function of order at least n in d. 
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If one now considers equations (15.1), (15.2), (15.3), (15.4), (15.5), (15.6), 
(15.7), (15.8), and (15.11), then the determinant of this sub-system of 
equations is —5184p*(p*+-1)(p?—1). Since p is a constant and p + 0, 
p # +1, this determinant is not zero, so that from these equations, one 
can explicitly and uniquely calculate ag, a, a9,..., xz. By Cramer’s Rule, 
these values have the form: 








OMPp . ss B, 
=r . . td - . 
Oo 
% - Eas ae |. Wars OO oe yk ee 
5184p%(p*+-1)(p?—1) 5184p%(p?+ 1)(p2—1) 
from which it follows that lim X; QO. for each i 0. l, e .. 8. This. 
d-0 


however, contradicts the assumption that an allowable 6th order analogue 
existed. Hence, a contradiction has been reached and the theorem is 
proved. 
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SUMMARY 


Generalizations to a rectangular grid are given of some fundamental lemmas and 
theorems of 8. Gerschgorin. The popular five-point technique for the numerical 
resolution of boundary-value problems associated with the Laplace equation is 
extended. The Dirichlet problem is replaced by the problem of solving a system 
of linear, algebraic equations, which, it is shown, has a unique solution. It is proved 
finally that the numerical solution converges to the analytic solution as the mesh 
constants converge to zero. 


1. Introduction 

In 1930 S. Gerschgorin proved some fundamental theorems germane to 
the numerical solution of elliptic differential equations (1). Since then, 
the vast advances in the field of computational machinery have added 
to the importance of that paper. One present need, and the intent here, 
is to generalize part of Gerschgorin’s work, which was accomplished for 
the square grid, to the rectangular grid. 

Let G, then, be a closed, bounded, simply connected plane region whose 
interior is denoted by R and whose boundary curve is denoted by 8. 
Let g(x,y) be defined and continuous on S. The Dirichlet problem is to 
produce a function u(z,y) such that both (a) and (6), as follow, are 
satisfied: 

(a) on G, u(x, y) satisfies the Laplace equation, that is: 

U,,+U,, = 9, (1.1) 

(b) on S, u(x, y) = g(x,y). 

Under general conditions, there exists a unique solution (2 and 3), and 
only such problems will be considered here. However, the analytical 
determination of u(x, y) is quite another story from that of its existence 
and usually offers what are at present insurmountable problems. Hence 
the approach here is through numerical analysis. 


2. General method 


Let A and d be fixed positive constants and let (x), y,) be an arbitrary, 
but fixed, point of G. Denote by G, the set of all points of the form 
(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 1, 1959) 
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(%)+mh, y.+nd), contained in G, where m and n are integers. Two points 
(2,,y¥,) and (2, y,) of G, are called adjacent if and only if 


(a) % = %, Yo—-%1| = 4; 
or (6) Y%=Y, La—2X,| = h. 


The interior of G,, denoted by R,, is the set of all points of G, which 
have four adjacent points which also belong to G,. The boundary of G,, 
denoted by S,,, and called the lattice boundary, is defined by 


ee ae 


The numerical method, then, is as follows. Suppose G, consists of n 
points. Number these points in a one to one fashion with the integers 


1, 2, 3,..., m. Denote the coordinates of the point numbered k by (2,, y;,) 
and the unknown function u at (x,,¥,) by 


u(t, Y,) =u, for k= 1, 2, 3,..., n. 


Let (x;,y;) be an arbitrary point of S,. Approximate u; by g(x’, y’), 
where (x’,y’) is the nearest point of S to (2;,y;). If (x’,y’) is not unique, 
choose any one of the set of nearest points and use it. The problem of 
finding numerical approximations to u(z,y) on the lattice boundary is, 
though crudely done, adequate for present purposes. We require then 
that at each point (x;,y;) of R,, the function wu satisfy 


—2u(x;,y;)-4 : palmar thy) +u(a;,—h, y;)]+ 


oe ’ i -[u(x;, y;+d)+ U(x; Y¥i— d)| =x @, (2.1) 
-+-p* 


where h=pd (h>0,d>0,p> 0). (2.2) 


Application of (2.1) to each point of R, results in a system of linear 
equations which, when solved, yields the remaining numerical approxima- 
tions. What remains to be shown then is the derivation of (2.1), the fact 
that the linear system just described has a unique solution and general 


conditions under which the numerical solution converges to the analytic 
solution. 


3. Derivation of the difference analogue (2.1) 


The grid sizes, by the construction of G,, are A in the z-direction and 
d in the y-direction. The following notation is used throughout: 


Uy = U(z,y), Uy = ulr+h,y), uu, = u(x,y+d), 
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\" 
(ar, ytd}® 2 
(xc-h,y) OL lac+h,y) 
3 ey) 1 


(nya 








Fic. 1 
Us = u(x—h,y), Uy = u(xz,y—d) (see diagram), (3.1) 
omy 
Ann = [1/(m!n!)] aay (3.2) 
Since w is harmonic, it follows from (4, p. 426), that 

(4) Ago = —Agg (i) Asg=—2Ayy | 
(b) Ago = —$A,2 (j) Ags = —10Ag,5 
(c) Ag = —3Ag; (k) Ago _ —Aogg 
(d) Ago = Ag () As, = Ais (3.3) 
(e) As, = —Ajs (m) Ag, = 15Agy, ; 
(f) Ase, = —6Ag, (n) As,;= —¥PA,, 
(9) Aso = tArs (0) Aga = —15Ag. 
(h) Ay, = 5Ags } 





For future reference, we note by use of Taylor’s theorem and by (2.2): 
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Uy = Ago 
Uy = Apo tA Pd+Agop*d®+ Ago p*d?+ Agoptd*+ 
+A, p*d°+ Ag, p*d®+ O(a’) 
Ut, = Ago t+Ao, d+ Ao2d*?+Aggd?+ Aggd*+ Ags d°+ Ao, gd°+O(d") \. (3,4) 
ts = Agg—Ajopd+ Azo p*d?—A,, p*d®+ Ag, p*d*—A,, p*d?+ 
T Agy p*d®+-O(d") 
Uy = Aog—Ao, d+ Ap. d?—Ayyd*+ Ay, d*— Ag, d°+ Ang d®+ O(d') 











4 
Now let L(u) = > a, u,. (3.5) 
0 


Substitution of (3.4) into (3.5), and then substitution of (3.3) into this 
result, yield, after recombination: 


L(u) = Aoo(49+4,+-a,+a,+4,)+d[A, 9 p(a,—ay)+ Ay ,(4.—4,)]- 
+d*Ag of (4.+-44)—p*(a, +45) ]+4*[Ao.3(42—44)— (1/3)p* A, 2(4, —45)] 4 
+d*Ao 4 p*(a, +43) + (4, +4,4)]+-4*[Ao5(42—44) + (1/5) p° Ay (4, —4)]4 
-+-d*A, .{ (a,+-a,)— p*(a, +a;)|+ O(d’). (3.6) 

In order to make (3.6) of as high an order in d as possible, consider 

Ay +A,+A,+dg+a, = 0, 
a,—a, = U, 
Aag—a, = V, 


A y+, p*(a, +s) 0. 


If one lets a, 2, then the solution of this system is 
a, = a, = 1/(1+-p?), a, = a, = p*/(1+p?). 
Hence, (3.5) becomes 
a ] ¥ me 
L(u) 2tlp + 5 (Uy +Us) 4 5 (Ma + M4) (3.7) 
+p l+-p 
while (3.6) becomes L(u) O(d*). (3.8) 
Then, from (3.7) and (3.8), 
1 2 
By +5 (+g) ++ P’ = (Ug+U4) = O(d'). (3.9) 
+ p* 1+-p* 


Discarding the terms O(d*) in identity (3.9) and recalling the definitions 
Of Ug, Uy, Ug, Ug, Uy yield the difference equation approximation (2.1). 

Note, of course, that for a square grid, that is, for p = 1, (2.1) is 
equivalent to the usual five-point analogue of Laplace’s equation, that is: 


~4u,+ Uj, +Usgt+Ust+u, = 0. 
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4. Theorems on the numerical method described in section 2 

In this section let u(z,y) be the analytic solution of the Dirichlet 
problem under consideration and let U(z,y) be the solution of the 
numerical method described in section 2. This latter assumption means 
that if (x, y) is any ies of R,, then 


—2U,+- (U,.+U;)+> (U,+U,) = (4.1) 


? = 


where 
» = U(a,y), U, = Ula+h,y), U, = U(x, y+d), 
U; = U(x—h,y), U, = U(x, y—d). 

THEoreM 1. The solution of the system of linear equations which results 
by application of the process of section 2 is unique. 

Proof. It is sufficient to show that the determinant of the system of 
linear equations is not zero and this is done by demonstrating that the 
only solution of the homogeneous system which results by considering 
g(x,y) = 0 on S is the zero vector. Suppose then there exists a non- 
trivial solution for the homogeneous system. For some point of R,, 
U +0. Suppose, without loss of generality, VU > 0. Let the largest value 
M occur at (29, ¥). Then: 

2U, 


ea 


zp +0) +70 (Ut Uy) (4.2) 


and M=U,2>U, (= ae (4.3) 
Case 1. Suppose U, = U, = U, = U,. Then from (4.2) and (4.3), 
o= U=M. 
Case 2. Suppose not all of U,, U,, U3, Uy are equal. Then there exists 
a maximum. Suppose that U, is the maximum. Hence 

U,=U,-t, U,s=U,—-t, U,= U,—t, (4.4) 


l 


where ¢,, fs, t, are non-negative and at least one is positive. Substitution 
of (4.4) into (4.2) yields 

ou, — 2U, — . tp? t,p® 

Wier t1+p? 1+p? 147? 
from which it follows that U, > Up, by the nature of f,, t;, t,. Hence, (4.3) 
is contradicted. A similar contradiction is reached if U,, U;, or Uy is 
selected as the maximum. Hence, in any case, U, = U, = M. 

Applying the same argument, a finite number of times, one may con- 
tinue this result and show that U’ at a boundary point is equal to M. But 
this is also a contradiction, since g(x,y) = 0 on S implies U = 0 at each 
point of S,. Hence the theorem is proved. 
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DEFINITION. Let 


h?+-d?2 


Liv(z,y)) = 8 | —20(2,4) + [v(ar+-h, y)+o(a—h, y)]+ 


1+) 
p* ) \ 
—___ [v(2, y+-d)-+ v(x, y—d)]}. 
Lemma |. If L{v] < 0 on R, and v > 0 on 8,, then v > 0 on R,. The 
proof follows easily as in (5, p. 156). 
Lemma 2. If —|L{v,]| > L[v,] on R, and }|v,| < v, on S,, then |v,| < v2 
on R,. The proof again follows easily as in (5, p. 156). 
Lemma 3. If |L[v]| < A on R, and if |\v| < BonS, and if r is the radius 
of any circle which contains G, then \v| < }(Ar*)+-B on R,,. 


Proof. Let 





ute, 9) | 40) Coe +a), 
| eee 

where (x—a)*+-(y—b)? = r? is the equation of a circle which contains G. 

Direct calculation yields: L(w) = —A. Also, w > B on S,. Now since 

L{v|| < A on R,, by assumption, and L[w] = —A, it follows that, on 

R,, |L{v]| < —L{w]. Also, |v} < B on S, and w > B on S,, so that 

w > |v! on S,. Hence —/|L[{v]| > L{w] on R, and w > jv} on S,. By 

Lemma 2, then, w > |v on R,, or |v} < w < }(Ar*)+B, on R,. Hence 
the lemma is proved. 


THEOREM 2. Jf u denotes the solution of the Dirichlet problem described 
in section 1, U denotes the solution of the linear system of the method described 
in section 2, then if (x;,y,;) is any point of R,, it follows that 





2(p24 2 

U(x;, y,)—u(ax;, y;)| - i d ) 4 Mh +d) (4.5) 
where 
M,= max [\t6(%,y)|, |Uo4(x,y))|], M, = max [|u,9(z, Y), Upr(%, y) |), 


(z,y)eG (z.yyeG 
r = radius of any fixed circle which contains G, 
om ny 
tan = 
ox™coy”™ 
Proof. Let Q = Liu|]—(u,,+4u,,). Substitution of the finite Taylor 
series expansions for u,, U, U3, uy in L[u] yields 


Q — WaolEr YAP , UaolE Wh® , taal Md? , tol, N2)d* 
4! , 4! : 4! : tees 
_ My(?+d2) 


ei 








Hence Q 











ON THE NUMERICAL SOLUTION OF DIRICHLET PROBLEMS 123 
Note, of course, that 
| L[u]+(u,.+u w)| = : | L{u)| — |Q|. 
Also, for - point of S,, U was taken as the value g(z’,y’) at the nearest 


point (x’,y’) on the boundary S, and g(z’,y’) = u(x’,y’) on S. Thereby, 
for any point (x,y) of S, we have 


U(x, y)—u(x, y)| = \g(z’,y’)—u(x, y)| = |u(z’, y')—u(z, y)| 
where (2’, y’) is a point of S and 
(x—2')?+(y—y')? < 2h, or (e@—2')P+(y—y')? < 2d. 
Therefore 


(U(x,y)—u(x,y)| = |u(x’, y’)—u(z, y’)+u(z, y’)—ulz,y) 


< |u,(Es, y’)(x—2’)|+ |u, (x, ns(y—y’)| 
< M,\x—2'|+|y—y'|] < My(h+d). (4.6) 
It must also be noted that L[U] = 0, by (4.1). Hence 
L{U0—u)| = |L[U]—L[u)| = |Lfu}| = |Q'. (4.7) 
Applying Lemma 3 to (4.6) and (4.7), then, on R, 
Ua) < ry TO) + agit d) = TO | asa) 


(4.8) 
which proves the theorem. 

THEOREM 3. Under the conditions of Theorem 2, the numerical solution 
U converges to the analytical solution u, as (h?+-d*) converges to zero. 


The proof for points of R, follows immediately from inequality (4.8), 
while the proof for points of S, follows directly from the numerical pro- 
cedure described for such points. 
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SUMMARY 
The motion of a pendulum with a periodic disturbance is examined. When higher 
powers of the displacement are retained there may be one or three periodic solutions 
according to the amplitude and period of the disturbing force. The elementary 
solution far from resonance corresponds to two different branches of the general 
solution, one of which ceases to exist when the resonance is sufficiently close. 


1. WHEN a small horizontal periodic force is applied to a pendulum, the 
displacement satisfies a differential equation of the form, with ¢« small, 
¢+o* sina = ecosxsint. (1) 
If x is neglected this has the elementary solution 
esint . 
‘a (2) 
and this is certainly a good approximation in most cases. However, when 
o* passes through | the coefficient of sint changes sign; this would suggest 
that the amplitude either passes through 0 or becomes infinite, and neither 
alternative is plausible—in fact the former does not satisfy the exact 
differential equation. 

The problem has been studied by E. W. Brown (1, 2, 3) who shows that 
when o? = | and higher powers of x are retained there is a periodic 
solution of (1) with amplitude of order e', and there is a suggestion that 
the amplitude for all o is at most of this order, but this has not been 
proved. In any case little is known about the behaviour in the range of 
transition. This paper concerns the existence and stability of periodic 
solutions. 

Case 1. o? >1. In the first place we notice that if o? > 1 and « = 0, 
there is an amplitude a that makes the period 27. In this motion z is an 
odd function of t+-«, where « is constant. Write it as X(t+-«). This is 
a first approximation to x; but if 


72+ ¢7(1—cos2), (3) 





P= 3 
Ej, = ecosxsint # dt, * 
0 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 1, 1959} 
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and at least for amplitude < 47 this increases indefinitely with ¢ unless 2, 
considered as a Fourier series, contains no term in sint and therefore x 
no term in cost. For a suitable first approximation to a periodic motion 
we must therefore take a = 0 or 7. For higher accuracy put 

x= X(t)+y (5) 
with y small. Then to the first order in « and y 

ij+o* cos X.y = ecos X sint. (6) 

The complementary function is of the form 

Au(t)+ Bvi(t), (7) 
where w(t), v(t) are an even and an odd function of t, and their Wronskian 
ué—iuv can be made 1. Then to order e 


t 
y = Au(t)+ Bo(t)— { {u(t)v(r)—u(r)v(t)}e cos X(r)sin + dr. (8) 
0 
The integral is an odd function of t; then we find that y(7) = y(—7), 
in) = y(—n), if 


A=0, Bv(m)— j {u(m)v(r)—u(r)v(7)}ecos X(z)sinr dr = 0. (9) 
0 


Then y is an odd periodic function. v(t) is not periodic, and v(7) + 0; 
hence B is determined and there are two periodic solutions of the form 

x = X(t)+c6(t), x= —X(t)+6(t) (10) 
where 6(¢) is an odd periodic function. 

Hence for o? > 1 there are three periodic solutions of (1), the elementary 
one and two large ones given by (10). For o? < 1 there is no analogue of 
the second pair, since any free motion would have a speed < o. 

If o? is near 1 the method fails because the amplitude of X is small 


and so is v(7). Hence B is the ratio of two small quantities and might not 
be small. 


Case 2. o? near 1. To investigate the transition we put 
oe—l=yp 
with » small, and neglect powers of x above the third. Then (1) becomes 
+a = fo%x—pxr+e(1—}2%)sint. (11) 
We shall put a=ez; p=elvy; d= 7x. (12) 
Then to the lowest powers of « 

+z = n(4z*—vz+8In 2), (13) 

where » is small but v is unrestricted. Substitute 
z = asint+csin 3t (14) 
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and assume c = o(a). The terms in sint and sin 3¢ give 
4a*®—}a*c—va+1 = 0, (15) 
8c = n(—3,a°+-4a*c—vr), (16) 


whence c is small of order na*; and we have to consider the cubic 


a’—8va+8 = 0. (17) 
v= 0 (o? = 1) gives a = —2, and the coefficient of sint in x is —2et, 


agreeing with Brown’s result. There is always a negative root. 
For v > 0 there are two positive roots if 


vy > 3.2¢ = 3.2-4, 


otherwise none. The critical value makes them coalesce at a = 2'. When 
v becomes large the smaller corresponds to the elementary solution, the 
larger to one of the large solutions. The negative root gives the other 
large solution. 

For v < 0 there are never three solutions. But with v << 0 the solution 
approximates to the elementary one a = 1/v. 

To investigate the variation of a with v we have 


da sa 


— — 18 
dy 3a*— 8p (18) 


For v < 0, a < 0 this is always negative. Hence as v increases from —<o 
to 0 the magnitude of a increases steadily to 2. As v proceeds to positive 
values the increase continues; 3a?—8v never vanishes for the negative 
root of (17), since the left side would be 


-Mya+8 > 8 


for v > 0, a < 0. Hence the elementary solution for v < 0 passes con- 
tinuously, with increasing amplitude, into the large negative solution for 
v > 0. For v < 3.2-' (nearly 1) this is the only periodic solution possible. 

For larger v there are two solutions with positive a, the smaller of which 
approximates to the elementary solution, but as v decreases through the 
critical value these disappear. If the motion has been periodic and near 
one of them, and » decreases, the further motion cannot be periodic at 
all; it may be described as an oscillation about the solution with negative 
a, but with period only approximately 27. The latter solution has 
a = —2.2§ = —3°3 nearly. 

Search has been made for other small solutions whose leading terms are 
of the form asin pt with p ~ 1; all attempts led to contradictions. 
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2. Stability of solutions 


If z is slightly disturbed from a periodic solution the disturbance may 
or may not ultimately become large. We write 


z = asint+{; (19) 

then {+0 = n(ja*sin%—v)f = n(}a*—v—}a* cos 2)f, (20) 
which is a Mathieu equation. It is known that the solutions of 

¢+(R—2q cos 2t)f = 0 (21) 

contain real exponential factors if R (near 1) lies between the values a,, b, 


a 


Unstable_ naa 
- —_ 


aa 
“ 
4 
/ 
"hss. sigue 


Stable 








Tic 





Fic. 1. Sketch of variation of a with v for the periodic solutions. 


that make an even and an odd periodic solution exist, namely (4, pp. 16, 


ati a, = 1+9—4q°..., (22) 

b, = 1—q—}¢".... (23) 

Here q=}a*, R= 1—}a*+», (24) 

a, = 1+}a*..., 6, = 1—}a’.... 

Then R = b, if ja* = v, and R = a, if ja? = v. The instability condition 
is ja® < » < ja’. 

From (17), v—}a* has the sign of a. Hence the condition is never satis- 

fied for vy positive or negative if a is negative. The solution with negative 


a is always stable. For positive a the first inequality is satisfied and 
stability depends on the sign of ja*—v. Since the vanishing of this 
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quantity is the condition for equal roots, when two positive values of a 
exist the smaller gives a stable solution and the larger an unstable one. 

The conclusions are therefore rather peculiar. If we imagine a sequence 
of systems with varying v, and start with large negative v (o® <1) the 
elementary solution starts by being a good approximation. As v increases 
the amplitude increases in magnitude and is 2e' when v = 0, but it con- 
tinues to increase as v becomes positive and may become large—but not 
at the values that make the elementary solution large. This corresponds 
(with a slight change of parameters) to the development of a large 
amplitude in a swing. 

If v is large and positive, and the solution is the large one with negative 
a, decreasing v will simply retrace the above process. But more usually 
the solution in this case will initially correspond to the small elementary 
one. With decreasing v the amplitude increases, but near v = 1 the solu- 
tion ceases to exist. With further decrease of v the motion will be a non- 
periodic one, oscillating slowly about the solution with negative a. 
Presumably damping will remove this oscillation in time. 

A term in 2? in the differential equation would not affect the main 
features of the solution. For on separation into harmonic parts it will give 
a constant and a term in cos 2t, neither of which would affect the coefficient 
of sint, and (17) would be unaltered. The solution would contain a constant 
and a term in cos 2¢, and would thus cease to be exactly an odd function 
of t, but the conclusions about the ranges of existence of the solutions and 
their stability would be hardly affected. 


In writing this paper I have benefited greatly from discussions with 
Professor J. E. Littlewood and Dr. Mary Cartwright. In particular they 
provided topological proofs of the existence of the solutions that are given 
approximately above. 
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